Final  Technical  Report 
AFOSR  Grant  No.  F49620-01-1-0074 

Mathematical  and  Computational  Framework  for 
Virtual  Fabrication  Environment  for 
Aircraft  Components 

Bama  Szabo1,  Sebastian  Nervi2  and  Daniel  Muntges3 
Center  for  Computational  Mechanics 
Washington  University,  Campus  Box  1129 
one  Brookings  Drive 
St.  Louis,  MO  63130-4899 


Acknowledgements 

This  work  was  sponsored  by  the  Air  Force  Office  of  Scientific  Research,  Air  Force  Materiel 
Command,  USAF,  under  grant  number  F49620-01-1-0074.  Valuable  experimental  data  were 
provided  by  the  Metallic  Processes  and  Prototyping  Laboratory  of  Boeing  Phantom  Works, 
St.  Louis,  Alcoa  Technical  Center,  Los  Alamos  National  Laboratory  and  Oak  Ridge  National 
Laboratory.  X-ray  diffraction  measurements  were  performed  at  WPAFB. 

Dr.  David  M.  Bowden,  Technical  Fellow,  Advanced  Manufacturing  Research  &  Develop¬ 
ment,  Mr.  Keith  Young,  Engineer/Scientist,  Metallic  Processes  and  Prototyping  Laboratory, 
Boeing  Phantom  Works,  St.  Louis,  Dr.  Mark  Newborn,  Alcoa  Technical  Center,  Dr.  Michael 
B.  Prime,  Technical  Staff  Member,  Los  Alamos  National  Laboratory  provided  advice  and 
assistance  in  this  investigation. 

The  U.S.  Government  is  authorized  to  reproduce  and  distribute  reprints  for  Governmental 
purposes  notwithstanding  any  copyright  notation  thereon. 

V 

Disclaimer 

The  views  and  conclusions  contained  herein  are  those  of  the  authors  and  should  not  be 
interpreted  as  necessarily  representing  the  official  policies  or  endorsements,  either  expressed 
or  implied,  of  the  Air  Force  Office  of  Scientific  Research  or  the  U.S.  Government. 

1The  Albert  P.  and  Blanche  Y.  Greensfelder  Professor  of  Mechanics. 

2Graduate  Research  Assistant. 

3Graduate  Research  Assistant. 


20050520  013 


2 


Contents 


1  Summary 

1.1  Air  Force  relevance  . 

1.2  Objectives . 

1.3  Accomplishments  .  . 

1.4  Archived  Publications 


4 

4 

7 

7 

9 


2  Verification  and  validation  11 

2.1  State  of  the  art  models .  11 

2.2  Working  models .  11 

2.3  Numerical  solution . _ . . .  12 

2.4  Reliability  of  computed  information . 12 

2.5  Why  is  numerical  accuracy  important? .  13 

2.6  Aspects  of  implementation . .  .  14 

3  Residual  stresses  16 

3.1  Path  independence .  16 

3.2  Validation  experiments . 18 

3.3  Sources  of  error  and  uncertainty .  19 

3.4  Results . 20 


4  Models  for  intersecting  shells  22 

4.1  Experiments . 23 

4.2  Formulation . 25 

4.3  Finite  element  spaces  . . . .  .  .  . .  29 

4.4  Working  models . 32 

4.5  Numerical  results  .  . .  35 

4.6  Summary  and  conclusions . .  .  40 

5  Generalized  plane  strain  42 


3 


1  Summary 


This  project  addressed  some  fundamental  questions  that  pertain  to  the  utility  of  the  methods 
of  numerical  mathematics  in  engineering  decision-making  processes.  The  basic  methodology 
investigated  concerns  the  problem  of  selection,  calibration  and  validation  of  mathematical 
models  and  thus  transcends  particular  applications.  The  principal  objective  was  to  perform 
an  investigation  on  the  development  of  cost-effective  methods  for  the  design  and  analysis  of 
unitized  metallic  airframe  components  that  account  for  the  effects  of  residual  stresses  and 
provide  sufficiently  reliable  information  to  serve  as  a  basis  for  certification  or  rejection  of 
manufactured  components. 


1.1  Air  Force  relevance 

This  project  addressed  some  of  the  key  requirements  identified  in  the  report  Uninhabited  Air 
Vehicles:  Enabling  Science  for  Military  Systems  (2000)  National  Materials  Advisory  Board, 
Aeronautics  and  Space  Engineering  Board  [1].  The  following  statements  are  quoted  from 
the  referenced  report: 


Analytical  tools 

“Research  should  be  initiated  to  integrate  design  and  analysis  models  and  methods  into  a  ver¬ 
satile  engineering  tool.  Current  analytical  models ,  such  as  the  finite  element  codes  developed 
by  the  National  Aeronautics  and  Space  Administration  (NASA )  and  structural  evaluation 
codes  developed  by  industry,  will  have  to  be  modified  and  improved  for  structural  design  and 
analysis  in  a  production  environment.  The  development  of  analytical  tools  should  also  include 
the  development  of  procedures  for  assessing  the  accuracy  and  reliability  of  model  predictions.” 
(p.  52) 


Characterization  and  Testing 

“Fundamental  research  should  be  undertaken  to  establish  potential  failure  modes  and  perfor¬ 
mance  levels  for  materials  and  structures  to  support  probabilistic  analysis  methods.  Current 
approaches  for  testing  materials  are  based  on  deterministic  design  methods  and  rely  on  ex¬ 
tensive  testing  at  the  subelement,  element,  subcomponent,  component,  and  full-scale  levels, 
using  a  building  block  approach  (NRC,  1996).  The  development  of  basic  property  relation¬ 
ships  and  potential  failure  modes  are  needed  for  implementing  probabilistic  design  approaches 
and  reducing  the  amount  of  large-scale  verification  testing  required.  ”  (p.  53) 


Simulation  Methods 

“Techniques  and  software  codes  should  be  developed  for  computational  simulations  of  struc¬ 
tural  responses  to  operational  environments  throughout  the  structures  lifetime  at  both  the 
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material  and  structural  levels.  Effective  analytical  simulations  would  enable  designers  to 
model  design  alternatives  without  developing  and  testing  expensive  prototypes,  resulting  in 
potentially  significant  reductions  in  developmental  costs.  ”  (p.  53) 

Design  Criteria 

“ Fundamental  research  on  critical  failure  modes  and  property  relationships  to  establish  mean¬ 
ingful  design  criteria  for  probabilistic  methods  should  be  undertaken.  Criteria  could  be  es¬ 
tablished  based  on  the  results  of  studies  on  the  relationship  between  the  conventional  safety 
factor  and  the  probabilistic  reliability  of  a  structure,  along  with  an  in-depth  survey  of  existing 
structures.”  (p.  53) 


Recommendation 

“To  support  the  development  and  introduction  of  probabilistic  methods  for  UAVs,  the  U.S. 
Air  Force  should  sponsor  research  on  (1)  analytical  tools,  (2)  characterization  and  testing, 
(3)  simulation  methods,  and  (4)  design  criteria.”  (p.  53) 

While  reference  [1]  addressed  design  needs  associated  with  the  development  of  UAVs,  a 
software  product  that  will  meet  its  recommendations  will  have  a  much  broader  applicability 
in  both  the  military  and  civilian  sectors:  The  value  of  information  acquired  through  various 
applications  of  the  techniques  of  numerical  mathematics  is  strongly  correlated  with  its  re¬ 
liability.  Reliable  information  will  serve  to  reduce  product  cycle  costs,  whereas  engineering 
decisions  based  on  unreliable  information  are  more  than  likely  to  cause  costly  problems.  Un¬ 
reliable  information  has  a  negative  economic  value:  Having  misleading  information  is  worse 
than  knowing  that  some  information  is  unavailable. 

It  is  well  known  that  finding  problems  late  in  a  product  cycle  is  generally  very  expensive 
in  terms  of  costs  and  performance.  The  failure  of  the  C-17  wing  in  a  static  test  in  October 
1992  is  a  case  in  point  [2].  The  costs  of  repair,  retrofitting  and  weight  penalties  typically 
can  be  traced  to  erroneous  decisions  made  in  the  design  phase. 

The  key  tQ  improved  reliability  in  numerical  simulation  is  systematic  application  of  two 
complementary  processes  called  verification  and  validation  (V&V).  Verification  is  concerned 
with  proper  approximation  of  the  exact  solution  of  a  mathematical  model,  called  working 
model.  The  goal  of  verification  is  to  ensure  that  the  exact  solution  is  approximated  with 
sufficient  precision,  so  that  decisions  based  on  data  computed  from  the  approximate  solution 
are  substantially  the  same  as  decisions  based  on  data  that  would  be  computed  from  the  exact 
solution.  Validation  is  concerned  with  proper  selection  of  the  working  model.  The  goal  of 
validation  is  to  ensure  that  simplifying  assumptions  incorporated  in  the  working  model  do 
not  affect  the  data  of  interest  significantly.  A  more  detailed  discussion  on  V&V  is  presented 
in  Section  2. 

With  very  few  exceptions,  V&V  processes  are  not  used  in  current  practice  and  hence  the 
full  potential  of  computer-aided  engineering  (CAE)  technology  is  not  being  realized.  The 
consequences  are  erroneous  and  sub-optimal  engineering  decisions,  unnecessary  and  poorly 
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conceived  physical  experimentation,  long  product  cycles  and  costly  maintenance.  This  very 
important  aspect  of  numerical  mathematics  is  receiving  an  increasing  amount  of  attention 
in  computational  fluid  dynamics  (CFD)  (see,  for  example,  [3],  [4],  [6],  [5])  but  has  not  yet 
received  much  attention  in  solid  mechanics. 

This  project  illustrates  the  benefits  that  can  be  derived  from  the  application  of  V&V  pro¬ 
cesses  to  thin-walled  unitized  structural  components  manufactured  by  means  of  high-speed 
machining  (HSM)  techniques.  A  test  article,  typical  in  complexity  of  unitized  structural 
components,  is  shown  in  Fig.  1. 


Figure  1:  Test  article  for  gridlock  experiments.  Some  of  the  panels  are  buckled  by  residual  stresses.  Repro¬ 
duced  with  permission  from  The  Boeing  Company. 

There  are  strong  economic  incentives  for  increasing  material  utilization  factors,  reducing 
machining  times  and  increasing  the  success  rate  for  machining  unitized  structural  compo¬ 
nents,  such  as  spars,  keel  beams,  bulkheads4.  A  reliable  virtual  fabrication  environment  is 
essential  for  controlling  the  costs  associated  with  fabrication  of  complex  aircraft  components 
and  improving  the  quality  of  the  design  process. 

There  is  a  need  for  specialized  software  products  that  provide  in-depth  simulation  capa¬ 
bilities  for  specific  classes  of  problems.  This  project  exemplifies  a  new  paradigm  in  numerical 
simulation  and  with  reference  to  an  important  problem  class:  the  design  and  certification  of 
unitized  metallic  airframe  components. 

4  Economic  estimates  have  been  developed  by  Dr.  David  M.  Bowden,  Technical  Fellow,  The  Boeing  Company,  St.  Louis, 
MO  63166  Tel.  314-232-1859,  david.m.bowden@boeing.com. 
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1.2  Objectives 


The  objectives  of  this  project  were  in  two  categories:  General  objectives  and  specific,  goal- 
oriented  objectives: 


1.2.1  General  objectives 

One  of  the  fundamental  questions  of  numerical  mathematics  is  how  mathematical  problems 
should  be  formulated  so  that  they  can  serve  as  reliable  and  accurate  representations  of 
physical  reality.  A  well-constructed  mathematical  model  should  provide  predictions  of  events 
or  states  of  systems  that  can  be  confirmed  consistently  by  physical  observations. 

The  general  objective  was  to  develop  guidelines  for  the  development  of  working  models 
for  the  simulation  of  the  structural  and  strength  responses  of  unitized  aircraft  components 
and  interpretation  of  physical  experiments  based  on  V&V  procedures. 

1.2.2  Specific  objectives 

The  specific  objectives  were  to  investigate  the  application  of  V&V  procedures  to  the  so¬ 
lution  of  specific  engineering  problems  of  substantial  significance  in  aerospace  engineering. 
This  requires  access  to  carefully  developed  experimental  data,  Based  on  the  availability  of 
experimental  data,  two  areas  were  selected  for  detailed  investigation: 

1.  Development  of  the  mathematical  and  computational  aspects  of  the  knowledge  base 
needed  for  the  creation  of  a  virtual  fabrication  environment  for  light  weight  aircraft 
and  spacecraft  components  manufactured  from  7050-T7451  aluminum  plate  stock  that 
will  make  it  possible  to  plan  the  fabrication  processes  so  that  the  incidence  of  re-working 
and  scrapping  of  partially  or  fully  manufactured  parts  is  substantially  reduced. 

2.  Investigation  of  working  models  for  the  analysis  of  thin-walled  structural  components 
fabricated  by  mechanical  milling. 

>- 

1.3  Accomplishments 

The  main  accomplishments  are  summarized  in  relation  to  the  general  and  specific  objectives 
in  the  following. 


1.3.1  General  objectives 

The  formulation  of  rules  for  design  and  fabrication  generally  involves  evaluation  of  experi¬ 
mental  data  with  the  objective  to  test  whether  a  working  model  or  hypothesis  concerning 
material  properties,  failure  criteria  or  boundary  conditions,  should  be  accepted  or  rejected. 
There  are  three  sources  of  error:  (a)  errors  in  the  working  model  or  hypothesis  being  tested, 
(b)  errors  in  the  numerical  approximation  and  (c)  errors  in  the  experiment.  Unless  the  errors 
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in  the  experiment  and  the  numerical  approximation  are  controlled,  it  will  not  be  possible  to 
judge  whether  the  working  model  or  hypothesis  should  be  rejected  or  not. 

1.  Destructive  physical  experiments,  performed  with  the  objective  to  determine  residual 
stresses,  involve  measurements  of  strains  and  displacements  in  specific  locations.  The 
distribution  of  the  residual  stress  is  inferred  from  the  observations  by  numerical  means, 
using  a  working  model  based  on  the  assumption  that  the  residual  stress  state  depends 
only  on  the  initial  stress  distribution  and  the  current  configuration  of  the  body.  This 
assumption  had  not  been  proven,  however.  A  proof  was  developed  under  this  project. 

2.  Working  models  used  for  the  purpose  of  interpretation  of  layer  removal  experiments 
are  typically  based  on  the  assumption  of  plane  strain  conditions.  In  reality,  however, 
the  conditions  under  which  these  experiments  are  performed  are  more  closely  approx¬ 
imated  by  generalized  plane  strain  than  plane  strain  conditions.  The  errors  in  the 
working  model  associated  with  the  use  of  plane  strain  rather  than  generalized  plane 
strain  conditions  were  estimated  by  theoretical  and  numerical  means. 

3.  The  choice  of  a  working  model  is  based  on  expert  opinion  formulated  with  the  aid  of 
virtual  and  physical  experimentation.  It  was  shown  that  in  many  cases  the  use  of  virtual 
experimentation  has  advantages  over  physical  experimentation. 

1.3.2  Specific  objectives 

1.  Physical  experiments  were  performed  and  independently  obtained  experimental  data, 
provided  by  Los  Alamos  National  Laboratory,  were  analyzed  in  an  effort  to  determine 
residual  stresses  induced  by  the  manufacturing  process  of  2  to  4-inch  7050-T7451  alu¬ 
minum  plate  stock.  The  experimental  arrangement  is  shown  in  Fig.  2. 

2.  The  distortion  of  thin-walled  test  articles  typical  of  unitized  aircraft  components  made 
of  7050-T7451  aluminum  was  predicted.  A  working  model  incorporating  the  assump¬ 
tions  of  the  linear  theory  of  elasticity  was  used  for  this  purpose.  The  results  indicated 
that  the  magnitude  of  distortion  that  can  be  attributed  to  the  residual  stresses  caused  by 
the  manufacturing  process  of  7050-T7451  plate  stock  is  not  substantial.  Consequently, 
either  the  working  model  is  not  suitable  for  the  intended  purpose,  or  the  large  distor¬ 
tions  observed  in  manufacturing  are  caused  by  machining-induced  stresses.  In  order  to 
determine  the  reason  for  the  differences,  validation  experiments  had  to  be  performed. 

3.  Validation  experiments  were  performed  with  the  objective  to  identify  the  effects  of  the 
bulk  stresses  (the  stresses  caused  by  the  manufacturing  process  of  the  plate).  The  results 
were  consistent  with  ’blind’  predictions  based  on  the  working  model.  Since  machining- 
induced  stresses  can  be  removed  by  chemical  milling,  the  results  indicate  that  high¬ 
speed  machining  should  be  used  for  the  fabrication  of  unitized  aircraft  components 
from  7050-T7451  aluminum  plate  stock  to  slightly  larger  then  the  intended  dimensions, 
then  chemical  milling  should  be  used  to  obtain  the  final  dimensions. 
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4.  Experimental  data  obtained  for  intersecting  thin- walled  cylindrical  shells  was  analyzed, 
utilizing  a  hierarchic  modeling  concept.  It  was  concluded  that  in  the  case  of  thin- 
walled  components,  subjected  to  loads  in  the  linearly  elastic  range,  complicated  physical 
experiments  are  not  warranted.  Virtual  experimentation  is  sufficient  for  the  purposes 
of  strength  analysis  of  such  components. 


(a) 


Figure  2:  Experiments  performed  at  Washington  University,  (a)  The  milling  machine  and  controls  (Cincin¬ 
nati  Milarcon  Sabre  machining  Centre),  (b)  Instrumented  specimens,  (c)  The  milling  tool:  3/4  inch  end 
mill  (two  flutes)  corner  radius:  0.120  inches  (3.05  mm). 


1.4  Archived  Publications 

The  following ‘papers  and  theses  that  document  various  aspects  of  this  effort  have  either 
been  published  or  are  in  the  process  of  being  published: 

Szabo  B  and  Actis  R.  On  the  importance  and  uses  of  feedback  information  in  FEA,  Applied 
Numerical  Mathematics.  52  (2005)  219-234. 

Szabo  B,  Duster  A  and  Rank  E.  The  p-version  of  the  finite  element  method.  In:  E.  Stein, 
R.  de  Borst,  T.J.R.  Hughes  (Eds.)  Encyclopedia  of  Computational  Mechanics.  Vol.  1, 
Fundamentals,  Ch.  5.  John  Wiley  &  Sons,  London  2004. 

Babuska  I  and  Szabo  B.  On  the  generalized  plane  strain  problem  in  thermoelasticity.  -  Sub¬ 
mitted  to  Comput.  Methods  Appl.  Mech.  Engng.  Under  review. 

Szabo  B  and  Muntges  D.  Procedures  for  the  verification  and  validation  of  working  models 
for  structural  shells.  Journal  of  Applied  Mechanics.  Paper  No.  JAM-04-1223.  In  press. 
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An  extended  version  of  this  paper  is  available  in  the  Proceedings,  2004  ASME  International 
Mechanical  Engineering  Congress,  Anaheim,  CA  Paper.  No.  IMECE2004-60014  November 
2004. 

Liebschutz  RI.  Simulation  of  residual  stresses  in  aluminum  alloy  using  p- version  finite  element 
methods.  MS  Thesis,  Washington  University,  June  2002. 

Muntges  D.  Verification  and  validation  of  working  models  for  thin  cylindrical  shells.  MS  the¬ 
sis,  The  Henry  Edwin  Sever  Graduate  School,  Washington  University,  St.  Louis,  December 
2004. 

Nervi  S.  A  mathematical  model  for  the  estimation  of  the  effects  of  residual  stresses  in  alu¬ 
minum  plates.  D.Sc.  Thesis,  Washington  University,  May  2005. 
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2  Verification  and  validation 


The  statement:  “  Effective  analytical  simulations  would  enable  designers  to  model  design 
alternatives  without  developing  and  testing  expensive  prototypes,  resulting  in  potentially  sig¬ 
nificant  reductions  in  developmental  costs.”  in  reference  [1]  emphasizes  the  desirability  of 
replacement  of  physical  testing  programs  with  numerical  simulation.  This  is  possible  if  the 
numerical  simulation  methods  are  sufficiently  reliable.  In  this  section  the  problem  of  quality 
assurance  in  numerical  simulation  is  discussed. 

The  main  elements  of  numerical  simulation  and  the  associated  errors  are  illustrated 
schematically  in  Fig.  3.  The  goal  of  CAE  is  to  obtain  reliable  quantitative  information 
concerning  some  physical  system.  To  achieve  this  goal,  a  working  model  has  to  be  defined 
and  a  numerical  solution  obtained.  These  points  are  discussed  in  the  following  sections. 


r  Errors  < 
JnterprotatlonJ 


to?  \ 
itatlon/ 


Figure  3:  The  main  elements  of  numerical  simulation  and  the  associated  errors. 


2.1  State  of  the  art  models 

The  fullest  possible  mathematical  representation  of  a  physical  system  or  object,  utilizing 
state  of  the  art  methods  and  information,  is  called  a  state  of  the  art  model.  State  of 
the  art  models  are  generally  very  complicated  and  the  information  sought  usually  can  be 
acquired  with*  simpler  models.  For  these  reasons  simplified  models,  called  working  models, 
are  formulated. 

2.2  Working  models 

Working  models  are  mathematical  problems  formulated  with  the  objective  to  capture  the 
essential  characteristics  of  state  of  the  art  models  with  the  expectation  that  engineering 
conclusions  based  on  them  will  be  substantially  the  same  as  if  the  state  of  the  art  model  had 
been  used. 

Associated  with  each  working  model  is  a  modeling  error  (indicated  schematically  in 
Fig.  3).  At  present  the  selection  of  a  working  model,  and  hence  the  control  of  modeling 
errors,  is  largely  left  to  the  experience  and  judgment  of  analysts  and  designers.  However, 
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as  the  complexity  of  the  systems  increases,  reliance  on  experience  and  judgment  becomes 
increasingly  impractical.  Expert  knowledge  must  be  aided  by  virtual  and  physical  exper¬ 
imentation.  Selection  of  the  appropriate  working  model  for  a  given  problem  is  the  most 
important  decision  in  any  CAE  application. 

Evaluation  and  modification  of  a  working  model  through  virtual  and  physical  experimen¬ 
tation  and  correlation  of  computed  information  with  experimental  data  is  called  validation. 
The  goal  of  validation  is  to  ensure  that  the  model  is  a  sufficiently  accurate  mathematical 
description  of  the  physical  system  or  process  it  is  supposed  to  represent  [3],  [4],  [6].  Valida¬ 
tion  involves  calibration  and  prediction.  The  determination  of  physical  properties  and  other 
model  parameters  through  correlation  with  experiments  is  called  calibration.  The  use  of  a 
model  to  foretell  the  state  of  a  physical  system  under  conditions  for  which  the  model  has 
not  been  calibrated  is  called  prediction. 

It  is  important  to  note  that  successful  prediction  of  the  outcome  of  one  or  more  physical 
experiments  does  not  prove  that  the  working  model  will  provide  correct  predictions  for 
other  experiments  or  applications.  It  shows  only  that  the  working  model  meets  necessary 
conditions  for  acceptance.  The  reliability  of  working  models  can  be  established  to  a  high 
degree  of  certainty  only  through  carefully  designed  validation  experiments  [4],  [6]. 


2.3  Numerical  solution 

The  solution  of  a  working  model  is  approximated  by  numerical  means,  most  frequently  by  the 
finite  element  method  (FEM)  [7].  It  is  necessary  to  ensure  that  the  data  of  interest  computed 
by  numerical  means  are  within  acceptable  error  bounds  with  respect  to  their  counterparts 
corresponding  to  the  exact  solution  of  the  working  model. 

Determination  of  the  accuracy  of  data  computed  from  the  approximate  solution  of  a 
particular  mathematical  model  is  called  verification.  In  a  verification  process  accuracy  is 
understood  to  be  with  respect  to  the  exact  solution  of  the  working  model,  not  with  respect 
to  physical  reality.  In  the  terminology  introduced  in  Fig.  3,  verification  is  concerned  with 
estimation  and  control  of  the  errors  of  approximation  and  determination  that  the  program 
is  functioningvas  intended  by  the  developer. 


2.4  Reliability  of  computed  information 

Information  generated  by  numerical  means  will  be  reliable  only  if  the  state-of-the-art  errors, 
the  modeling  errors  and  the  errors  of  approximation  are  sufficiently  small  and  no  conceptual 
errors,  or  errors  of  interpretation,  or  programming  errors  have  occurred.  Proper  selection 
of  a  working  model  depends  on  the  goals  of  computation.  Commonly  occurring  conceptual 
errors  include  attempting  to  compute  data  that  are  inconsistent  with  the  formulation  of  the 
working  model.  The  detection  of  programming  errors  is  usually  quite  difficult,  necessitating 
careful  checking  of  the  program  for  the  problem  class  and  the  data  of  interest. 

The  state  of  the  art  errors,  the  modeling  errors  and  the  errors  of  approximation  are 
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benign  errors  in  the  sense  that  these  are  controllable  errors.  State-of-the-art  errors  can  be 
reduced  by  experimentation;  modeling  errors  and  errors  of  discretization  can  be  reduced  by 
automated  or  user-guided  adaptive  processes.  Conceptual  errors,  errors  of  interpretation 
and  programming  errors  are  malignant,  in  the  sense  that  they  render  the  results  of  analysis 
useless,  or  worse,  misleading. 


2.5  Why  is  numerical  accuracy  important? 

It  is  often  asked:  “If  we  do  not  know  the  loading  data  with  sufficient  accuracy  then  why  do  we 
need  to  be  concerned  with  the  accuracy  of  the  numerical  solution ?”  In  answering  this  question 
we  consider  two  important  areas:  The  application  of  design  rules  and  formulation  of  design 
rules.  It  is  shown  in  the  following  that  neither  proper  formulation  nor  proper  application  of 
design  rules  is  possible  without  estimation  and  control  of  the  numerical  accuracy. 


2.5.1  Formulation  of  design  rules 

The  formulation  of  design  rules  involves  evaluation  of  experimental  data  with  the  objective  to 
test  whether  a  working  model  or  hypothesis  concerning  material  properties  or  failure  criteria 
should  be  accepted  or  rejected.  There  are  three  sources  of  error:  (a)  errors  in  the  working 
model  or  hypothesis  being  tested,  (b)  errors  in  the  numerical  approximation  and  (c)  errors 
in  the  experiment.  Unless  the  errors  in  the  experiment  and  the  numerical  approximation  are 
controlled,  it  will  not  be  possible  to  judge  whether  the  working  model  or  hypothesis  should 
be  rejected  or  not. 


2.5.2  Application  of  design  rules 

Design  and  design  certification  involve  application  of  existing  design  rules,  established  by 
various  codes  and  conventions  of  professional  practice.  The  design  conditions  and  criteria 
to  be  used  are  mandated  by  the  cognizant  authorities.  There  is  no  uncertainty  concerning 
the  loading  conditions  or  other  types  of  excitation  that  must  be  considered,  or  the  physical 
properties  to  be  used.  The  design  rules  are  typically  stated  in  the  form  of  required  minimum 
factors  of  safety  (FS): 

FS:=*-*r  T  (!) 

^maxv^-EX/ 

where  $iim  >  0  is  the  limiting  (not  to  exceed)  value  of  some  functional  (such  as  maximum 
stress)  and  ^max(uEx)  >  0  is  the  exact  value  of  the  same  functional  corresponding  to  the 
exact  solution  of  the  mathematical  model.  The  designer’s  responsibility  is  to  ensure  that  the 
applicable  design  rules  are  followed. 

We  will  denote  by  $raax(«F£)  the  value  of  $max  computed  from  the  finite  element  solution. 
Let  us  suppose  that,  owing  to  numerical  errors,  it  is  possible  to  guarantee  only  that  the 
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relative  error  is  less  than  r: 


l^max^Ex)  ~  $max(uFE)| 


<  T  0  <  T  <  1 


<f>ma X(UEX)  -  ' 

i.e.,  <&max(uFE)  may  underestimate  $ma x(uex)  by  lOOr  percent.  Therefore  we  have: 

^max  (uex)  <  7—^-  \  ^max  (ufe) 

(1  ~T) 

on  substituting  this  expression  into  eq.  (1),  we  have: 


FS  <  (1  -  t) 


$ma x(ufe) 


Denoting  the  factor  of  safety  to  be  used  in  conjunction  with  the  data  generated  by  finite 
element  analysis  by  FSfe • 

FS"  ■=  «) 


we  have  from  eq.  (2): 


FSfe  > 


1  —  r 


This  result  shows  that  to  compensate  for  numerical  errors  in  the  computation  of  <hma x(ufe), 
it  is  necessary  to  increase  the  required  factor  of  safety  to  FS/(  1  —  r).  For  example,  if  the 
accuracy  of  $ma x(«fe)  can  be  guaranteed  to  20%  (i.e.,  r  =  0.20)  then  the  factor  of  safety 
must  be  increased  by  25%  to  compensate  for  numerical  errors.  Since  the  factor  of  safety  was 
already  chosen  conservatively  to  account  for  various  uncertainties,  the  economic  penalties 
associated  with  using  an  increased  factor  of  safety  generally  far  outweigh  the  costs  associated 
with  guaranteeing  the  accuracy  of  the  data  of  interest  to  within  a  small  relative  error  (say 
5%). 

It  is  necessary  to  specify  the  acceptable  error  tolerance  in  numerical  analysis  and  to  verify 
that  the  error  is  not  larger  than  the  specified  tolerance  even  when  that  tolerance  is  large. 


2.6  Aspects  of  implementation 

In  the  early  papers  and  books  on  finite  element  analysis  it  was  customary  to  write  strain  in 
the  following  form 

{e}  =  [D][N]{a}  =  [J3]{a} 

where  [N]  is  the  matrix  of  element-level  basis  functions  and  {a}  are  the  coefficients  of  the 
basis  functions  (see,  for  example  [10],  [11])-  Therefore  the  element-level  stiffness  matrix  was 
written  as: 

[Ke}=  [  [B]T[E][B]dV 
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rather  than  in  the  form: 


[*.]=  /  (P][iV])T[B][J5][W]<iV. 

-  Jiic 

In  this  way  the  formulation  of  the  working  model,  represented  by  the  operator  [D],  the  poly¬ 
nomial  degree  of  the  elements  and  the  mapping  functions  which  are  related  to  discretization, 
were  combined.  This  lead  to  the  development  of  large  element  libraries  with  separate  element 
types  for  each  working  model  and  polynomial  degree,  such  as  the  constant  strain  triangle,  the 
Argyris  triangle,  Bell’s  triangle,  Hsieh-Clough-Tocher  triangle  for  plate  bending  (a  descrip¬ 
tion  and  analysis  of  which  is  available  in  [13]),  various  shallow  shell  elements,  etc.  There  are 
elements  in  commercial  computer  codes  for  which  the  details  of  the  formulation  are  not  made 
available.  The  infrastructures  of  these  codes  were  designed  to  support  control  of  discretiza¬ 
tion  errors  through  mesh  refinement,  using  elements  of  fixed  and  low  polynomial  degree,  but 
do  not  provide  for  systematic  control  of  modeling  errors.  For  a  discussion  on  shell  elements 
we  refer  to  [12]. 

In  the  hierarchical  view  the  formulation  and  the  discretization  must  be  separated  in  order 
to  provide  a  suitable  computational  framework  for  selecting  a  suitable  working  model  and 
for  verifying  the  computed  information.  Furthermore,  working  models  must  be  defined  so 
that  they  have  consistent  meaning  when  passing  from  one  model  to  another. 

Considering  control  of  discretization  errors,  once  again  a  hierarchic  framework  is  needed. 
The  goal  of  computation  is  to  approximate  some  functionals  $i(uEx)>  i  —  1, 2, . . . ,  to  within 
a  specified  tolerance.  Of  course,  the  exact  solution  is  generally  not  known,  however  certain 
properties  of  the  exact  solution  are  known  a  priori.  It  is  necessary  to  construct  a  sequence 
of  finite  element  spaces  Si(Q)  C  C  •  •  •  C  Sk(Q)  C  •  •  *  and  compute  the  corresponding 

finite  element  solutions  uEE  and  from  the  finite  element  solutions  the  functionals 
k  —  1,2,...  This  has  very  substantial  impact  on  the  design  of  an  implementation.  Conven¬ 
tional  implementations  were  not  designed  in  this  way.  For  this  reason  they  are  not  suited 
well  for  systematic  control  of  the  errors  of  idealization  or  the  errors  of  discretization. 
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3  Residual  stresses 


In  this  section  we  address  the  specific  objective  of  developing  a  model  that  can  predict,  with 
reasonably  accuracy,  the  distortion  of  a  part  machined  from  large  metal  plates,  in  particular 
from  50.8  mm  (2  in)  to  101.6  mm  (4  in)  thick  7050-T7451  aluminum  plates,  caused  by  residual 
stress  introduced  in  the  manufacturing  process  of  the  plates.  The  pattern  of  residual  stress 
in  this  aluminum  plate  is  complicated  because  following  quenching,  the  plates  are  stretched 
in  order  to  reduce  the  magnitude  of  residual  stresses.  Since  the  data  that  characterize  the 
yield  behavior  of  the  material  vary  through  the  thickness,  a  complicated  pattern  of  residual 
stress  develops,  which  is  illustrated  in  Fig.  5. 


Figure  4:  Residual  stress  distribution  in  the  transverse  direction. 

The  determination  of  residual  stresses  is  by  destructive  experiments.  In  these  experiments 
displacements  and/or  strains  are  measured  from  which  the  residual  stress  state  is  inferred. 
The  experimental  arrangement  for  the  determination  of  residual  stresses  by  the  wide  slot 
method  is  shown  in  Fig.  2.  The  underlying  theory  is  based  on  the  assumption  that  the 
principle  of  superposition  holds,  from  which  the  property  of  path  independence,  discussed 
in  the  following  section,  can  be  proven. 


3.1  Path  independence 

Let  us  denote  the  initial  domain  of  a  body  by  Q0>  its  boundary  points  by  and  the  initial 
stress  state  by  cr?  .  The  initial  stress  state  satisfies  the  equations  of  equilibrium  cr?  ^  =  0 
on  Q0  and  the  stress-free  boundary  conditions,  i.e.  a^rij  =  0  on  5O0  where  rij  is  the  unit 
normal. 

Let  us  now  introduce  a  cut  rl5  that  produces  the  domain  Cl\,  and  a  second  cut  r2  to 
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produce  a  domain  ft2.  Denote  the  displacement  field  following  the  first  (resp.  the  second) 
cut  by  Uf  (resp.  Uf).  The  notation  is  shown  in  Fig.  5. 


Figure  5:  Path  independence.  Notation. 
Theorem  3.1  Uf  depends  on  crfj  and  Q2  but  not  on  Uf  or  Qj. 


Proof:  Uf  satisfies 


B(Uf,Vi )  =  —  [  afjTijVi  dS  V  Vi  e  E(Qi) 

JTi 


where  B(Uf,vf)  is  the  bilinear  form  corresponding  to  the  principle  of  virtual  work,  Vi  is  a 
test  function  and  E(Qi)  is  the  energy  space. 

Denote  by  Sfj  the  stress  field  corresponding  to  Uf  on  fii  and  let  uf  be  the  displacement 
field  on  02  corresponding  to  the  tractions  Sf^rij  acting  on  F2,  therefore 

B(uf,  Vi)  =  f  SfjTijVi  dS  V  Vi  €  E( 02).  (5) 

Jr2 

The  displacement  field  corresponding  to  the  second  cut  is  Uf  —  uf  +  Uf,  where 

B(Uf,  Vi)  =  -  f  (a°.  +  Sfjnj)vi  dS  V  Vi  e  E(n2).  (6) 

JTi 

Adding  eq.  (5)  and  eq.  (6)  we  find  that 

B(Uf,  vf)  —  f  ofjUjVi  dS  V  t*  €  £^(^2)  (7) 

*■  *'r2 

which  proves  the  theorem. 


Remark  3.1  Here  we  have  considered  an  unconstrained  body,  hence  the  displacements  are 
understood  to  be  displacements  up  to  rigid  body  displacements.  Of  course,  a  part  cannot 
be  machined  without  constraints.  However,  if  the  effect  of  the  constraints  is  elastic,  and  the 
constraint  conditions  do  not  change  during  the  machining  process,  then  the  result  obtained 
in  eq.  (6)  remains  unchanged,  therefore  Theorem  3.1  holds  under  these  conditions  also. 


Remark  3.2  We  have  assumed  that  when  the  body  is  cut,  the  cutting  operation  will  not 
introduce  stresses.  However  in  mechanical  milling  additional  residual  stresses  are  introduced. 
The  magnitude  of  these  stresses  are  substantial  only  in  a  thin  layer  near  the  surface.  This 
problem  is  discussed  in  the  following  section. 
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3.2  Validation  experiments 


Validation  experiments  were  performed  at  the  Metallic  Processes  and  Prototyping  Labora¬ 
tory  of  Boeing  Phantom  Works,  St.  Louis.  The  objective  of  the  experiments  was  to  test 
whether  the  linear  working  model  meets  necessary  criteria  for  acceptance.  The  working 
model  is  based  on  the  assumption  that  if  the  stresses  introduced  by  the  machining  process 
were  negligibly  small  then,  in  accordance  with  Theorem  1,  the  configuration  of  a  workpiece 
after  machining  is  predictable  from  the  shape  of  the  part  and  the  residual  stress  distribution 
in  the  plate  from  which  it  was  machined. 

The  validation  experiments  were  planned  jointly  by  the  research  groups  at  the  Center 
for  Computational  Mechanics  of  Washington  University  and  Boeing  Phantom  Works  in  St. 
Louis.  A  detailed  description  of  the  experimental  procedure  is  available  in  [9].  It  was  decided 
to  machine  four  test  articles  having  a  Z-shaped  cross-section  from  a  single  4  inch  (101.6  mm) 
thick  7050-T7451  aluminum  plate.  The  nominal  dimensions  of  the  test  articles  are  shown 
in  Fig.  6.  The  thicknesses  of  the  test  articles  ranged  from  0.04  to  0.12  inches.  In  addition, 
test  specimens  were  cut  from  the  same  plate  and  the  residual  stresses  were  determined  by 
the  wide  slot  method  [8].  These  experiments  were  similar  to  those  performed  at  Washington 
University  earlier. 


NOMINAL  DIMENSIONS 


Dimension 

mm 

inches 

Thickness  ti 

3.048 

0.12 

Thickness  t2 

2.032 

0.08 

Thickness  £3 

1.524 

0.06 

Thickness  <4 

1.016 

0.04 

Flange  width  li 

40.64 

1.60 

Flange  width  li 

25.4 

1.00 

Web  w 

50.8 

2.00 

Fillet  radius  17 

3.048 

0.12 

Length  L 

914.4 

36.00 

v  Figure  6:  Cross-section  and  nominal  dimensions  of  the  test  articles. 

Some  machining-induced  stresses  are  invariably  present  following  mechanical  milling  op¬ 
erations.  Based  on  existing  experimental  evidence  it  is  expected  that  the  machining-induced 
stresses  are  confined  to  within  a  small  layer  at  the  surface  of  the  material  and  therefore  their 
effects  on  distortion  would  increase  as  the  thickness  of  the  test  articles  decreased. 

Machining-induced  stresses  can  be  eliminated  by  chemical  milling.  Assuming  that  the 
principle  of  superposition  holds,  and  the  data  obtained  for  the  aluminum  plate  are  reliable,  it 
should  be  possible  to  predict  the  deformed  configuration  of  a  test  article  once  the  machining- 
induced  stresses  were  removed.  Referring  to  Fig.  6,  the  two  specimens  of  thicknesses  £i  and 
t2  were  chemically  milled  to  thicknesses  i3  and  U  in  an  effort  to  eliminate  the  residual 
stresses  introduced  by  the  machining  process  and  to  have  a  direct  comparison  between  the 
configurations  of  test  articles  of  the  same  thickness,  one  of  which  is  affected  by  the  machining- 
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induced  stresses,  the  other  is  not.  The  four  Z-shaped  test  articles  attached  to  the  plate  by 
thin  tabs,  are  shown  in  Fig.  7. 


Figure  7:  The  four  Z-shaped  test  articles  after  machining.  The  test  articles  are  supported  by  tabs. 


The  working  model  proposed  for  predicting  the  distortion  caused  by  bulk  stresses  was 
solved  numerically  using  the  finite  element  method,  and  the  solution  was  verified  using  hierar¬ 
chic  sequences  of  finite  element  spaces.  The  model  was  tested  with  the  objective  to  determine 
whether  it  meets  necessary  conditions  for  acceptance.  This  test  involved  correlation  with 
experiments  using  a  ‘blind’  prediction:  The  analyst  (S.  Nervi)  made  the  prediction  without 
having  seen  the  outcome  of  the  experiments  and  the  person  conducting  the  experiments  (K. 
Young)  was  not  aware  of  the  predicted  configuration. 

The  configuration  of  the  Z-shaped  test  articles  was  measured  following  mechanical  and 
chemical  milling  on  three  faces  using  a  contact  probe  with  a  measurement  error  of  ±0.0127  mm 
(±0.0005  in).  The  measured  data  were  adjusted  to  ensure  that  the  same  coordinate  system 
was  used  in  measuring  the  configurations  prior  to  and  following  chemical  milling. 

3.3  Sources  of  error  and  uncertainty 

When  comparing  the  results  of  physical  experiments  with  data  obtained  by  numerical  sim¬ 
ulation,  it  is  necessary  to  consider  the  differences  between  the  real  physical  object  and  its 
mathematical  description. 
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3.3.1  Thickness  variation 


In  the  mathematical  model  uniform  thickness  was  assumed.  In  the  test  articles  there  were 
small  variations  in  thickness  that  lead  to  variations  in  the  displacement  field. 

3.3.2  Variation  of  residual  stress  over  the  plate 

The  wide  slot  method  used  for  obtaining  the  bulk  stresses  has  inherent  sources  of  error  and 
uncertainties,  such  as  the  effect  of  the  machining-induced  stresses  caused  by  the  cutting  tool 
used  for  creating  the  slot,  and  the  error  associated  with  the  measured  strains  which  will 
lead  to  errors  in  the  computed  value  of  the  bulk  stresses.  Furthermore,  the  bulk  stresses 
were  obtained  from  samples  taken  from  a  different  location  than  the  test  articles  and  it 
was  assumed  that  the  bulk  stresses  were  uniform  over  the  plate.  Since  we  considered  a 
linear  material  law,  variations  in  the  magnitude  of  the  residual  stresses  will  lead  to  similar 
variations  in  the  magnitude  of  the  predicted  distortion. 


3.3.3  Material  properties 

The  material  properties  used  were  standard  values,  no  experiments  were  conducted  to  deter¬ 
mine  the  modulus  of  elasticity  ( E )  and  Poisson’s  ration  (v).  The  method  used  for  computing 
the  residual  stresses  is  based  on  a  linear  material  law  for  which  the  computed  residual  stresses 
will  “adjust”  any  actual  difference  between  the  real  value  and  the  value  used  in  the  numerical 
simulation.  Therefore,  differences  between  the  actual  values  of  the  material  properties  and 
the  ones  used  in  the  simulations  will  not  affect  the  computed  distortion  but  they  will  affect 
the  computed  residual  stress  after  material  removal. 


3.3.4  Measurement  errors 

The  configuration  of  the  test  articles  was  measured  by  supporting  each  test  article  in  three 
points  using  spherical  contact  probes.  The  test  articles  were  bonded  to  the  contact  probes. 
The  coordinates  were  measured  using  an  automated  contact  probe,  which  approached  the 
test  article  along  the  normal  to  its  surface  until  contact  was  made.  Since  contact  is  necessary 
in  order  to  obtain  a  measurement,  the  probe  could  affect  the  displacement  when  the  test 
article  is  thin.  Furthermore,  the  effect  of  gravity,  not  considered  in  the  numerical  simulation, 
could  also  have  a  small  effect  on  the  configuration  of  the  test  articles. 


3.4  Results 

The  distortion  of  the  test  articles  increased  with  decreasing  thickness.  This  is  consistent 
with  the  expectation  that  the  effects  machining-induced  stress  are  strongly  correlated  with 
thickness. 
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The  predicted  and  measured  configurations  of  the  chemically  milled  test  articles  were  very 
similar  and  substantially  different  from  the  configuration  of  test  articles  that  were  not  milled 
chemically.  The  predicted  distortion  was  somewhat  smaller  than  the  observed  distortion. 
This  can  be  explained  by  the  differences  between  the  test  articles  and  their  mathematical 
representation,  the  variations  in  thickness  are  likely  to  be  the  main  reason  for  the  differences. 
Details  are  available  in  [8]. 
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4  Models  for  intersecting  shells 


In  structural  shells,  such  as  pressure  vessels  and  piping  systems,  there  are  regions  where  the 
kinematic  assumptions  incorporated  in  classical  shell  models  do  not  hold.  These  are  the 
neighborhoods  of  nozzles  and  shell  intersections,  support  attachments,  stiffeners,  cut-outs, 
and  regions  where  the  curvature  abruptly  changes.  Prom  the  point  of  view  of  engineering 
design  and  analysis,  these  are  the  typical  regions  of  primary  interest,  however. 

In  formulating  mathematical  models,  shell-like  structures  should  be  viewed  as  fully  three- 
dimensional  solid  bodies  that  allow  the  imposition  of  restrictions  on  the  transverse  variation 
of  displacement  vector  components  in  certain  regions.  The  imposition  of  such  restrictions  re¬ 
duces  a  three-dimensional  problem  to  a  two-dimensional  one  in  the  case  of  plates  and  shells; 
to  a  one-dimensional  problem  in  the  case  of  beams  and  arches.  An  important  practical  prob¬ 
lem  is  to  identify  the  proper  formulation  for  a  particular  application,  that  is,  a  formulation 
that  simplifies  the  problem  by  dimensional  reduction  without  affecting  the  data  of  interest 
significantly. 

In  our  terminology  a  working  model  is  a  mathematical  model  of  some  physical  system  or 
process  comprised  of  a  formulation  and  data  pertaining  to  geometric  description,  material 
properties,  loading  and  constraints.  Associated  with  a  working  model  is  an  exact  solution. 

Common  to  all  mathematical  models  based  on  the  principles  of  continuum  mechanics  are 
the  conservation  of  momentum  (in  static  problems  the  equations  of  equilibrium),  the  strain- 
displacement  relations  and  constitutive  laws.  A  particular  application  involves  selection  of 
a  working  model. 

Any  working  model  should  be  viewed  as  a  special  case  of  a  general  model,  called  state 
of  the  art  model.  By  definition,  a  state  of  the  art  model  accounts  for  all  physical  laws  that 
pertain  to  the  system  or  process  being  modeled.  Selection  of  a  particular  working  model 
involves  acceptance  of  certain  limitations  that  may  or  may  not  be  justified,  given  the  goals  of 
computation.  In  general  it  is  not  known  a  priori  whether  the  limitations  of  a  particular  work¬ 
ing  model  are  justified.  Therefore  it  is  necessary  to  perform  tests  a  posteriori  to  determine 
whether  a  particular  working  model  is  suitable  with  respect  to  the  goals  of  computation. 
If  the  tests  shpw  that  a  working  model  is  inadequate  then  a  more  comprehensive  working 
model  must  be  chosen  and  the  process  repeated. 

In  principle,  one  could  continue  this  process  indefinitely  until  all  relevant  attributes  of 
the  physical  system  or  process  being  modeled  are  taken  into  account.  In  reality  however 
the  more  comprehensive  a  working  model,  the  more  descriptive  data  that-  characterize  the 
geometric  description,  material  properties  and  boundary  conditions  are  needed.  Often  the 
needed  parameters  are  not  available,  or  would  be  difficult  to  determine.  The  cognitive 
as  well  as  the  stochastic  uncertainties  increase  with  the  complexity  of  the  working  model. 
Therefore  a  highly  detailed  working  model  is  generally  not  feasible.  Consequently  is  some 
degree  of  uncertainty  associated  with  any  attempt  to  model  a  physical  system  or  process  by 
mathematical  means.  On  the  other  hand,  much  useful  information  can  be  developed  through 
systematic  evaluation  of  sequences  of  idealized  working  models  with  reference  to  the  data  of 
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interest. 

Errors  associated  with  the  choice  of  model  are  called  errors  of  idealization  whereas  errors 
associated  with  finite  element  or  other  numerical  approximation  of  the  exact  solution  Of  the 
model  are  called  errors  of  discretization.  The  process  by  which  errors  of  idealization  are 
controlled  is  called  validation.  The  process  by  which  errors  of  discretization  are  controlled  is 
called  verification.  In  this  paper  control  of  both  the  errors  of  idealization  and  discretization 
are  discussed  with  reference  to  a  model  problem  for  which  detailed  experimental  information 
of  high  quality  is  available. 

The  model  problem  discussed  in  this  paper  is  the  first  of  four  models  investigated  by 
experimental  and  analytical  means  at  Oak  Ridge  National  Laboratory  in  the  1970’s  [15], 
[14].  The  test  articles  consisted  of  two  thin- walled  circular  cylindrical  shells  intersecting  at 
right  angles.  The  cylindrical  shells  were  made  of  carbon  steel.  The  test  articles  are  idealized 
representations  of  shell  intersections  that  occur  in  piping  and  pressure  vessel  systems  in  the 
sense  that  there  were  no  transitions,  fillets  or  reinforcements  at  the  junctions.  The  stated 
objective  of  the  experiments  was  to  obtain  stress  data  for  comparison  with  the  predictions 
of  thin  shell  theory. 

A  hierarchic  family  of  working  models  that  allow  virtually  arbitrary  representation  of  the 
transverse  variation  of  the  displacement  vector  components  is  described  in  the  following.  The 
lowest  member  of  the  hierarchy  is  similar  to  the  Naghdi  shell  model.  The  highest  member 
is  the  fully  three-dimensional  representation.  Selection  of  a  working  model  for  a  particular 
application  is  based  on  feedback  information  that  includes  the  results  of  tests  conducted 
with  the  objective  to  determine  whether  the  data  of  interest  are  dependent  on  the  restrictive 
assumptions  incorporated  in  the  model. 


4.1  Experiments 

Physical  experiments  were  performed  at  the  Oak  Ridge  National  Laboratory  in  the  1970’s 
[14],  [15].  The  goal  of  the  experiments  was  to  determine  whether  the  classical  model  for  shells, 
known  as  the  Novozhilov-Koiter  model,  discretized  by  an  assembly  of  flat  plate  elements,  is 
capable  of  predicting  strains  in  the  vicinity  of  the  intersection  of  two  cylindrical  shells.  Four 
carbon  steel  test  articles  were  manufactured  and  instrumented  with  great  care.  A  detailed 
analysis  of  the  first  test  article,  based  on  hierarchic  modeling  techniques,  is  available  in  [18]. 

The  first  test  article  was  made  by  welding  two  carbon  steel  pipes  then  carefully  machining 
the  weldment  to  the  test  dimensions.  The  test  article  was  annealed  at  several  points  in  the 
machining  process.  The  experimental  arrangement  is  shown  in  Fig.  8.  The  horizontal  part 
is  called  the  cylinder,  the  vertical  part  is  called  the  nozzle.  The  length  of  the  cylinder  was 
39.0  in  (991  mm).  The  length  of  the  nozzle,  measured  from  the  point  of  intersection  of  the 
centerline  of  the  nozzle  with  the  centerline  of  the  cylinder  was  19.5  in  (495  mm).  The  outside 
diameter  of  the  cylinder  (resp.  nozzle)  was  10  in  (254  mm)  (resp.  5.0  in  (127  mm)).  The 
intended  wall  thickness  of  the  cylinder  (resp.  nozzle)  was  0.1  in  (2.54  mm)  (resp.  0.05  in 
(1.27  mm)). 
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Figure  8:  Experimental  arrangement.  Reproduced  with  permission  from  Oak  Ridge  National  Laboratory. 
4.1.1  Constraint  conditions 

As  shown  in  Fig.  8,  the  right  end  of  the  cylinder  was  rigidly  clamped  to  a  heavy  flat  plate 
bolted  to  a  frame.  Small  flanges  were  machined  into  the  ends  of  the  cylinder  and  nozzle  to 
support  the  seal  and  the  clamping  forces.  Heavy  loading  fixtures  were  attached  on  the  left 
end  of  the  cylinder  and  the  end  of  the  nozzle  to  provide  seal  and  seating  for  the  application 
of  forces. 


4.1.2  Loading  conditions 

A  total  of  thifteen  (13)  load  cases  that  included  pressure  loading  and  axial  forces,  shear 
forces  and  moments  were  investigated.  The  forces  and  moments  were  applied  to  the  cylinder 
and  the  nozzle  through  hydraulic  rams  acting  through  load  cells.  The  pressure  loading 
was  applied  by  means  of  a  hydraulic  fluid.  In  order  to  compensate  for  the  weight  of  the 
hydraulic  fluid,  a  counterbalancing  force  was  applied  to  the  fixture  at  the  free  end  of  the 
cylinder  through  a  cable  that  is  visible  in  Fig.  8. 

For  all  13  load  cases  the  load  was  applied  in  increments  of  20  percent  of  the  full  load,  then 
decreased  to  zero  again  in  20  percent  decrements.  In  this  paper  only  one  of  the  loading  cases, 
the  pressure  load,  is  discussed.  Additional  load  cases  are  discussed  in  [18].  The  maximum 
value  of  the  pressure  was  50.0  psi  (344.8  kPa). 
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4.1.3  Measurements 


A  total  of  322  three-gauge  (Micro-Measurements  type  EA-06-030YB-120,  option  SE)  foil 
rosettes  were  bonded  on  the  inside  and  outside  surfaces  by  epoxy  adhesive  and  cured.  The 
gauges  in  the  rosettes  were  arranged  in  a  ‘Y’  pattern  (i.e.,  the  directions  of  measurement 
were  120  degrees  apart)  [14].  Details  of  the  instrumented  intersection  region  is  shown  in 
Fig.  9. 


Figure  9:  Detail  of  test  article  instrumented  with  strain  gauges.  Reproduced  with  permission  from  Oak 
Ridge  National  Laboratory. 

The  ratio  of  the  resistance  change  in  a  strain  gauge  to  the  Lagrangian  strain  causing  the 
change  is  called  the  gauge  factor.  The  gauge  factor  of  each  production  lot  is  determined  by 
sample  measurements  and  is  given  on  each  package  with  its  tolerance.  Typical  tolerances 
are  0.5  to  1.0  %.  No  tolerance  data  are  provided  in  references  [14],  [15]. 

4.2  Formulation 

The  formulation  of  working  models  for  structural  shells  is  a  large  and  complicated  subject 
that  cannot  be  discussed  in  sufficient  detail  here.  Only  a  brief  overview  of  the  main  points 
relevant  to  the  present  investigation  is  presented. 


4.2.1  Kinematic  assumptions 

A  structural  shell  is  characterized  by  a  surface,  called  mid-surface  a;,-,  and  the  thickness  t. 
Both  are  given  in  terms  of  two  parameters  ai,  <22: 

Xi  =  xt-(ai,0!2)>  t  =  t(ai,  a2). 

The  indices  for  a,-  take  on  the  values  i  =  1, 2  whereas  the  indices  of  x  range  from  1  to  3. 
Associated  with  each  point  of  the  mid-surface  are  three  basis  vectors.  Two  of  the  basis 
vectors  lie  in  the  tangent  plane: 

z,(i)  ._  i.(2)  9xi_ 

*  ‘  dax'  *  •  da2 

Note  that  b ^  and  are  not  necessarily  orthogonal.  The  third  basis  vector  bf^  is  the 
cross  product  of  b\^  and  bf\  therefore  it  is  normal  to  the  tangent  plane.  These  are  called 
curvilinear  basis  vectors.  The  normalized  curvilinear  basis  vectors  will  be  denoted  by  ea,ep, 
en.  The  Cartesian  unit  basis  vectors  corresponding  to  the  coordinates  Xi  will  be  denoted  by 
ex,  ey,  ez.  A  vector  u,  given  in  terms  of  the  curvilinear  basis  vectors,  denoted  by  U(„),  can 
be  transformed  to  Cartesian  coordinates,  denoted  by  U(x).  The  transformation  is 

u(x)  =  [fl]u(o)  (8) 

where  the  columns  of  the  transformation  matrix  [R]  are  the  unit  vectors  eQ,  ep,  en.  The 
displacement  vector  components  are  given  in  the  following  form: 

ma 

ua  :=  ^Ua\i(a,P)4>i{y) 
i= 0 
rap 

UP  :==  UP\i(a>  (9) 

»=0 
mn 

un  :=  ^2un]i(a, p)4>i(u) 

i= 0 
v 

where  u  is  the  independent  variable  in  the  direction  of  the  normal.  The  functions  ua|,-, 
u/9|tS  wn|i  are  called  field  functions,  the  functions  are  called  director  functions.  When 
the  material  is  isotropic  then  <&( u)  are  polynomials;  when  the  shell  is  laminated  then  4>i(u) 
are  piecewise  polynomials  (see,  for  example,  [22],  [19]).  Equation  (9)  represents  a  semi¬ 
discretization  in  the  sense  that  4>i(v)  are  fixed,  hence  the  number  of  dimensions  is  reduced 
from  three  to  two.  The  kinematic  assumptions  incorporated  in  a  particular  shell  model  are 
characterized  by  the  indices  (ma,  rrip,  mn).  The  lowest  member  of  the  hierarchy  is  the  model 
(1,  1,  0)  which,  from  the  point  of  view  of  kinematic  assumptions,  is  the  same  as  the  Naghdi 
shell  model  [20].  The  kinematic  assumptions  of  the  Novozhilov-Koiter  model  [21]  are  more 
restrictive  than  those  of  the  model  (1,  1,  0)  in  that  the  field  functions  uQ|x  and  up\\  are 
constrained  to  be  linear  combinations  of  the  first  derivatives  of  u„|o,  i.e.,  there  are  only  three 
independent  field  functions. 
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The  classical  development  of  shell  models  was  strongly  influenced  by  the  limitations  of 
the  methods  available  for  solving  the  resulting  systems  of  equations.  The  use  of  curvilin¬ 
ear  coordinates  allowed  the  treatment  of  shells  with  simple  geometric  description,  such  as 
cylindrical,  spherical  and  conical  shells  by  classical  methods,  subject  to  the  assumption  that 
the  thickness  of  the  shell  is  small  in  relation  to  its  other  dimensions.  Such  limitations  no 
longer  exist.  It  is  possible  to  formulate  the  problem  in  terms  of  either  the  curvilinear  or 
the  Cartesian  components  of  the  displacement  vector.  When  the  formulation  is  based  on 
the  curvilinear  (resp.  Cartesian)  components  of  the  displacement  vector  then  we  refer  to  the 
formulation  as  a  shell  (resp.  thin  solid)  formulation. 

In  the  following  we  will  be  concerned  with  the  thin  solid  formulation  only,  that  is,  the 
formulation  in  terms  of  the  Cartesian  components  of  the  displacement  field: 

m 

Ux  :=  y\x| i{a,  p)(j>i{v) 

»=o 

m  " 

uv  :=  53 P)Mv)  (10) 

i=0 

m 

uz  :=  53 
*=0 

Note  that,  in  the  case  of  thin  solid  models,  the  kinematic  assumptions  are  characterized  by 
the  single  index  m.  In  other  words,  the  transverse  variation  of  the  three  displacement  vector 
components  is  approximated  by  the  same  functions  (f>i{v),  i  =  0, 1, 2, . .  .m. 

Certain  advantages  and  disadvantages  are  associated  with  formulating  shell  models  in 
terms  of  the  Cartesian,  rather  than  the  curvilinear,  components  of  the  displacement  vector. 
The  advantages  are  that  continuity  with  other  bodies,  such  as  stiffeners,  are  easier  to  enforce 
and  implementation  is  simpler.  The  disadvantages  are  that  thin  solid  formulations  cannot  be 
applied  to  laminated  shells  unless  each  lamina  is  explicitly  modeled  or  homogenized  material 
properties  are  used;  the  number  of  field  functions  must  be  the  same  for  each  displacement 
component.  For  example,  the  (1,  1,  0)  shell  model  has  five  field  functions,  the  thin  solid 
model  characterized  by  m  =  1  has  six  field  functions. 

4.2.2  Linear  working  models  for  shells 

In  the  hierarchic  view  of  working  models  the  highest  (i.e.,  most  comprehensive)  model  ac¬ 
counts  for  all  possible  nonlinear  effects,  such  as  geometric  nonlinearities,  material  nonlinear¬ 
ities  and  mechanical  contact.  In  the  case  of  beams,  plates  and  shells  there  is  an  important 
subset  of  model  hierarchy  where  the  highest  model  is  the  fully  three-dimensional  problem  of 
linear  elasticity.  For  this  subset,  the  accuracy  of  the  exact  solution  of  a  working  model  is  un¬ 
derstood  to  be  in  relation  to  the  exact  solution  of  the  corresponding  fully  three-dimensional 
problem  of  linear  elasticity,  not  the  underlying  problem  of  solid  mechanics.  This  subset  is 
discussed  in  this  section. 
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The  vector  of  strain  tensor  components  corresponding  to  u  is  denoted  by  {e}.  The 
relationship  between  u  and  {e}  is  given  by  {e}  =  [D] u  where  [D]  the  differential  operator 
associated  with  small  strain  elasticity.  The  vector  of  stress  tensor  components  is  denoted  by 
{a}.  The  relationship  between  {a}  and  {e}  is  given  by  Hooke’s  law:  {a}  =  [i?]{e}  where 
[E\  is  the  elastic  material  stiffness  matrix.  Traction  vectors  acting  on  dQ.  are  denoted  by  T. 
The  virtual  work  of  internal  stresses  is  defined  as  follows: 

B(u,v)  :=  f  ([D]v)T[E][D]ndQ  (11) 

J  n 

and  the  virtual  work  of  external  forces  is  defined  by 

F(v)  :=  [  F  •  v dSl  +  [  T  -vdS+  [  ([D]v)t[E]{c}t  dCl  (12) 

Jtl  Jd(l  JCl 

where  {c*}  :=  {ct  Ct  ct  0  0  0}T,  c*  is  the  coefficient  of  thermal  expansion  and  r  =  t(x,  y,  z ) 
is  the  change  in  temperature  from  a  reference  temperature. 

The  energy  space  is  defined  by  E(Q)  :=  (u  |  B(n,  u)  <  C  <  oo}  where  C  is  some  positive 
constant.  The  energy  norm  is  defined  by 

\  1/2 

JB(u,u)J  .  (13) 

The  generic  form  of  the  principle  of  virtual  work  is  stated  as  follows:  Find  u  €  E(Q.)  such 
that 

B(u,v)  =  F(v)  for  all  v  e  E(Q).  (14) 

Specific  statements  of  the  principle  of  virtual  work  depend  on  the  boundary  conditions. 
For  details  we  refer  to  [7].  In  the  generic  case,  and  whenever  the  prescribed  displacement 
constraints  are  insufficient  to  prevent  all  possible  rigid  body  displacements,  the  solution  of 
eq.  (14)  is  unique  up  to  rigid  body  displacements. 

There  is  an  important  difference  between  the  classical  shell  models,  such  as  the  Novozhilov- 
Koiter  and  the  Naghdi  model,  and  the  hierarchic  shell  and  thin  solid  models.  In  using  hi¬ 
erarchic  modqls  the  goal  is  to  approximate  the  fully  three-dimensional  solution,  hence  the 
stress-strain  law  is  that  of  the  three-dimensional  theory  of  elasticity.  In  index  notation: 

Oij  =  A  Sijekk  +  2  y.ei:j  (15) 

where  aij,  Cij  are  the  Cartesian  stress  and  strain  tensors,  respectively,  A  and  /i  are  the  Lame 
parameters  and  is  the  Kronecker  delta.  Incorporated  in  the  stress-strain  relationship  of 
the  Naghdi  and  the  Novozhilov-Koiter  models  is  the  assumption  that  the  stress  component 
normal  to  the  mid-surface  is  zero,  a  condition  which  is  the  limiting  case  of  the  fully  three- 
dimensional  solution  with  respect  the  thickness  approaching  zero.  The  Naghdi  shell  model 
yields  the  correct  solution  in  the  limit  t  — >-  0,  however  neither  the  hierarchic  shell  model  (1, 
1,  0)  nor  the  thin  solid  model  m  —  1  does,  unless  Poisson’s  ratio  is  zero.  The  Naghdi  model 
is  not  a  member  of  the  hierarchic  sequence  of  models  but  rather  an  extension  of  the  sequence 
for  small  t  values.  Hierarchic  shell  (resp.  thin  solid)  models  characterized  by  mn  >  3  (resp. 
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m>  3)  give  the  correct  limit  solution  with  respect  to  t  -*  0  and  are  said  to  be  asymptotically 
consistent. 

Remark  4.1  In  the  case  of  shells  the  distinction  between  the  notions  of  mathematical  model 
and  its  discretization  is  blurred  by  conventions  in  terminology:  It  is  customary  to  refer  to 
the  various  shell  formulations  as  theories  or  models.  In  fact,  the  hierarchic  models  are  semi¬ 
discretizations  of  the  fully  three-dimensional  model.  Therefore  modeling  errors  that  can  be 
attributed  to  the  choice  of  indices  ( ma ,  mp,  mn )  in  the  case  of  hierarchic  shell  models,  and 
the  index  m  in  the  case  of  thin  solid  models,  are  related  to  discretization  rather  than  model 
definition. 


4.2.3  Nonlinear  working  models  for  shells 

The  model  hierarchy  must  account  for  nonlinear  effects.  This  large  and  important  topic  is 
not  within  the  scope  of  this  paper.  For  discussion  and  examples  we  refer  to  [16],  [17],  [23]. 


4.3  Finite  element  spaces 

The  accuracy  of  the  finite  element  solution  is  determined  by  the  finite  element  space.  Finite 
element  spaces  are  constructed  by  partitioning  the  solution  domain  Q  into  finite  elements. 
A  partition  will  be  denoted  by  A,  the  number  of  elements  of  the  partition  by  M(  A)  and  the 
fcth  element  by  Qk.  Typically  is  mapped  from  a  corresponding  standard  element  fist  by 
smooth  mapping  functions  Q^. 

Ideally,  finite  element  spaces  are  constructed  by  adaptive  methods  such  that  the  tolerance 
criteria  for  the  data  of  interest  are  satisfied.  The  initial  finite  element  mesh  should  be  laid  out 
utilizing  a  priori  information  concerning  the  regularity  of  the  exact  solution.  For  example, 
in  the  neighborhood  of  singular  points  and  lines  the  mesh  should  be  graded  in  geometric 
progression  with  a  fixed  common  factor  when  p-extension  is  used  [7].  If  h-extension  is  used 
then  radical  grading  is  optimal  [24].  In  the  case  of  plate  and  shell  models  the  presence  of 
boundary  layers  has  to  be  taken  into  account. 

V 

4.3.1  The  finite  element  space  used  in  the  ORNL  investigation 

At  the  time  of  the  ORNL  investigation  the  treatment  of  shell  models  by  the  finite  element 
method  was  in  its  very  early  stages  of  development.  Although  the  investigators  were  aware 
of  some  contemporary  work  on  curved  shell  elements,  shells  were  commonly  approximated 
by  flat  plate  elements  and  most  of  the  available  experience  was  with  those  elements.  For 
this  reason  the  investigators  decided  to  use  flat  plate  elements  for  the  purpose  of  analyzing 
the  shell  intersection  problem  [14].  The  Hsieh-Clough-Tocher  (HCT)  triangular  element  [25] 
was  chosen  for  the  approximation  of  the  displacement  component  normal  to  the  mid-surface 
of  the  shell.  The  constrained  linear  strain  triangle  (CLST)  was  used  for  approximating  the 
membrane  components.  A  brief  description  follows. 


29 


The  HCT  triangle  is  a  composite  element,  comprised  of  three  sub-triangles.  On  each 
sub-triangle  an  incomplete  cubic  polynomial  approximation,  comprised  of  9  terms,  is  used. 
The  polynomials  are  chosen  so  that  along  the  external  edges  the  normal  derivative  varies 
linearly.  Therefore  there  are  9  coefficients  per  sub-triangle.  C°  continuity  is  enforced  for 
the  sub-elements  by  constructing  basis  functions  corresponding  to  the  3  nodal  displacements 
and  6  rotations  for  each  sub-triangle,  as  indicated  in  Fig.  10(a)  where  the  circles  represent 
transverse  displacements  and  the  arrows  represent  first  derivatives  (rotations  in  the  sense  of 
the  arrows).  The  sub-triangles  are  assembled,  which  is  equivalent  to  satisfying  C°  continuity 
over  the  triangle. 


Figure  10:  (a)  The  12  degrees  of  freedom  HCT  triangle,  (b)  Composite  non-planar  quadrilateral  element 
assembled  from  four  HCT  triangles. 

At  this  point  there  are  3  internal  degrees  of  freedom,  indicated  in  Fig.  10(a)  by  the 
closed  circles  and  arrows,  and  9  external  degrees  of  freedom,  indicated  by  the  open  circles 
and  arrows.  In  order  to  satisfy  exact  and  minimal  C1  continuity,  the  continuity  of  the 
normal  derivatives  along  the  internal  edges  is  enforced  leading  to  3  constraint  equations  that 
establish  a  relationship  between  the  sets  of  internal  and  external  degrees  of  freedom.  Using 
these  constraint  equations  the  internal  degrees  of  freedom  are  eliminated.  Thus  the  HCT 
triangle  has  9  degrees  of  freedom:  three  displacements  and  three  rotations  in  each  coordinate 
direction. 

Non-planar  quadrilateral  plate  elements  assembled  from  four  HCT  elements  were  used  for 
approximating  the  displacement  vector  components  normal  to  the  shell  surface.  A  typical 
element  is  shown  in  Fig.  10(b).  Since  the  center  node  generally  does  not  lie  in  the  same 
plane  as  the  vertex  nodes,  there  is  a  third  rotation  component,  not  present  in  the  constituent 
triangles.  A  third  rotation  component  is  also  present  in  the  assembly  of  the  quadrilateral 
elements  into  the  stiffness  matrix,  since  adjacent  elements  are  generally  not  co-planar.  The 
usual  treatment  is  that  the  rotation  components  are  transformed  into  a  Cartesian  system,  the 
origin  of  which  lies  on  the  shell  surface  at  the  node  and  one  axis  is  coincident  with  the  normal 
to  the  surface.  The  rotation  component  in  the  direction  of  the  normal  is  usually  neglected. 
This  causes  various  problems,  for  details  we  refer  to  [14],  [15],  [26].  For  a  discussion  on  low 
order  shell  elements  we  refer  to  [12]. 

The  in-plane  (membrane)  components  of  the  displacement  vector  were  approximated 
by  a  similar  assembly  of  triangles.  These  vector  components  were  approximated  over  each 
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component  triangle  by  quadratic  polynomials  constrained  so  that  the  variation  is  linear 
over  the  external  edges.  This  is  known  as  the  constrained  linear  strain  triangle  (CLST).  A 
quadrilateral  membrane  element  has  two  degrees  of  freedom  at  each  of  its  five  nodes  [14]. 
The  finite  element  space  is  the  span  of  the  assembled  of  HCT  and  CLST  triangles  shown  in 
Fig.  11. 


END  PLATE 


Figure  11:  Finite  element  mesh  used  in  the  ORNL  investigation.  There  are  649  nodal  points  of  which  25 
nodal  points  lie  on  the  intersection.  Source:  ORNL-DWG  69-10664R.  Reprinted  with  permission  from  Oak 
Ridge  National  Laboratory. 


The  investigators  recognized  two  sources  of  error:  The  errors  caused  by  using  plate 
elements  rather  than  curved  shell  elements  and  by  using  a  limited  number  of  elements. 
Implied  is  the  tacit  assumption  that,  as  the  number  of  elements  is  increased,  the  approximate 
solution  will  converge  to  the  exact  solution  of  the  Novozhilov-Koiter  shell  model.  This  has  not 
been  proven.  There  is  no  guarantee  that  the  sequence  of  finite  element  solutions  obtained 
by  h-extension  will  converge  to  the  exact  solution  of  the  Novozhilov-Koiter  model  of  this 
problem. 


Remark  4.2  ‘The  definition  of  the  HCT  triangle  given  in  [13]  is  different  from  the  one 
described  here  in  that  all  sub-triangles  are  complete  polynomials  of  degree  3  in  which  case 
there  are  12  degrees  of  freedom,  the  nine  degrees  of  freedom  shown  in  Fig.  10(a)  plus  the 
first  derivative  in  the  direction  of  the  normal  at  the  mid-point  of  each  side. 


4.3.2  The  finite  element  spaces  used  in  the  present  investigation 

In  the  present  investigation  the  thin  solid  formulation  implemented  in  the  finite  element 
analysis  software  product  StressCheck5.  was  used.  The  finite  element  mesh,  consisting  of 
188  elements,  is  shown  in  Fig.  12.  The  nozzle  and  the  shell  were  partitioned  into  hexahedral 

5StressCheck  is  a  trademark  of  Engineering  Software  Research  and  Development,  Inc.,  St.  Louis,  Missouri. 
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elements,  the  intersection  region  was  partitioned  into  hexahedral  and  pentahedral  elements. 
The  hexahedral  elements  were  mapped  from  the  standard  hexahedron 

!={({;■?,  0 1  III  <x;  W<1,ICI<1}  (16) 

by  smooth  mapping  functions.  For  details  on  mapping  procedures  we  refer  to  [27]. 


Figure  12:  The  188-element  mesh  used  in  the  present  investigation. 

The  standard  polynomial  spaces  used  in  the  present  investigation  are  known  as  anisotropic 
trunk  spaces.  For  p>l,  l<q<p  these  spaces  are  defined  by: 

STfSl^J^spanK^'r.^r,  iffC, 

k,E  =  0, 1, 2, ... ,  k  +  E  <p,  m  =  0, 1, 2 . . . ,  q).  (17) 

For  definitions  of  the  isotropic  trunk  space  on  the  standard  hexahedron  <Sfr(fi ^),  and  the 
standard  pentahedron  we  refer  to  [7].  The  finite  element  space  S(Q)  is  the  span 

of  the  mapped  basis  functions  defined  on  the  standard  elements,  subject  to  the  exact  and 
minimal  continuity  requirement  of  the  formulation.  Referring  to  eq.  (10),  q  =  m,  corresponds 
to  the  thin  solid  formulation  characterized  by  the  index  m.  For  additional  details  we  refer 
to  [16],  [23],  [28]. 


4.4  Working  models 

The  quality  of  a  working  model  is  determined  by  the  proximity  of  its  exact  solution  to  the 
exact  solution  of  the  state  of  the  art  model,  more  precisely,  the  proximity  of  functionals  that 
would  be  computed  from  the  exact  solution  of  the  state  of  the  art  model  and  the  working 
model.  Here  we  examine  a  sequence  of  working  models  from  the  point  of  view  of  their  ability 
to  approximate  physical  observations  and  data  computed  from  physical  observations  through 
simple  transformation. 
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The  objective  of  the  ORNL  experiments  was  to  investigate  how  well  thin  shell  theory  is 
capable  of  predicting  the  stress  distribution  in  intersecting  cylindrical  shells  in  the  neighbor¬ 
hood  of  the  intersection.  The  investigators  had  in  mind  the  Novozhilov-Koiter  shell  model 
Only.  This  model  is  a  reasonable  choice  for  representing  the  smooth  parts  of  the  intersecting 
shells.  However,  given  that  the  data  of  interest  are  the  strain  values  in  the  vicinity  of  the 
intersection,  neither  the  kinematic  assumptions  nor  the  material  properties  (plane  stress) 
incorporated  in  this  model  are  valid  in  that  region.  Whereas  the  mathematical  problem  of 
the  intersecting  cylinders  for  the  Novozhilov-Koiter  shell  model  is  well  defined  [29],  it  is  not 
a  good  representation  of  the  physical  problem. 

The  working  models  employed  in  the  present  investigation  were  also  based  on  the  theory 
of  elasticity  but  differ  from  the  Novozhilov-Koiter  shell  model  in  kinematic  assumptions 
and  material  properties,  as  described  in  Section  4.2.  The  intersection  region  was  treated 
as  a  three-dimensional  elastic  region  for  each  model.  The  size  of  the  intersection  region 
(characterized  by  the  dimensions  ds  and  dn  shown  in  Fig.  13)  is  fixed  for  each  working 
model.  Elsewhere  the  kinematic  assumptions  of  the  thin  solid  formulation  were  used  with 
m  —  q—  1,2, ... 


Figure  13:  Mesh  detail  at  the  intersection  region  and  location  of  the  strain  gauges  in  the  plane  of  symmetry 
nearest  to  the  intersection.  188-element  mesh. 

The  problem  of  selecting  a  working  model  with  respect  to  the  goals  of  computation  of  the 
ORNL  experiments  is  understood  as  follows:  A  particular  working  model  based  on  the  thin 
solid  formulation  is  judged  to  meet  the  necessary  conditions  for  acceptance  if  the  equivalent 
(von  Mises)  stresses  computed  in  the  gauge  locations  do  not  differ  by  more  than  r  percent 
from  the  corresponding  equivalent  stresses  computed  from  the  fully  three-dimensional  model. 
The  choice  of  tolerance  is,  of  course,  arbitrary,  however  the  accuracy  required  of  the  numerical 
solution  depends  on  this  choice.  The  relative  error  in  the  numerical  solution  has  to  be  less 
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than  approximately  0.5r  percent. 

Using  the  188-element  mesh  shown  in  Fig.  12  and  the  anisotropic  spaces  Sffg(Q^)  with 
q  =  1,2, . . .,8  all  working  models  yielded  consistent  results  within  the  prescribed  tolerance 
of  r  =  2.5  %.  Therefore  the  simplest  model  ( q  =  1)  is  preferred. 

In  comparing  data  computed  from  working  models  with  physical  measurements,  it  is 
necessary  to  recognize  the  differences  between  the  mathematical  problem  being  solved  and 
the  physical  system  being  modeled.  These  differences  are  enumerated  in  relation  to  the 
ORNL  test  article  in  the  following. 

1.  Geometric  variations.  The  ORNL  investigators  were  careful  to  minimize  errors  in  man¬ 
ufacturing  the  test  article  however,  owing  to  unavoidable  machining  tolerances,  some 
variations  in  wall  thickness  and  the  other  dimensions  had  to  be  present.  Quoting  from 
reference  [14]; 

“A  careful  dimensional  inspection  of  the  machined  model  indicated  that,  despite 
the  care  taken  in  machining,  there  were  wall  thickness  variations  in  both  nozzle 
and  cylinder  with  the  nozzle  thickness  being  as  much  as  15  %  greater  (0.007  to 
0.008  in.  compared  with  the  nominal  0.050  in)  in  the  fourth  quadrant  than  in 
the  second.  ” 

In  the  working  models  it  is  assumed  that  the  shells  are  defined  by  perfect  cylindrical 
surfaces  and  constant  wall  thickness. 

The  intent  of  the  ORNL  investigators  was  to  manufacture  the  intersection  with  zero 
fillet  radius.  In  reality,  the  milling  tool  leaves  some  fillet,  see  Fig.  9.  In  the  working 
models  the  fillet  radius  is  zero. 

The  test  article  had  small  flanges  at  the  ends.  The  working  models  do  not  account  for 
those  flanges. 

2.  Variations  in  material  properties.  The  material  properties  assumed  by  the  ORNL  in¬ 
vestigators  are  the  nominal  elastic  constants  of  carbon  steel.  Modulus  of  elasticity: 
E  =  30  x  106  psi  (207  GPa);  Poisson’s  ratio:  u  =  0.3.  The  actual  elastic  constants  of 
the  test  Article  can  differ  from  the  nominal  values  by  a  few  percent.  There  are  no  data 
on  the  statistical  variations  of  the  modulus  of  elasticity  and  Poisson’s  ratio,  however  it 
is  reasonable  to  expect  that  the  mean  value  of  E  (resp.  v)  is  within  about  2  %  (resp.  5 
%)  of  the  nominal  value.  Therefore  systematic  as  well  as  random  errors  are  present  in 
making  comparisons  between  measured  strains  and  strains  computed  from  the  working 
models. 

The  relationship  between  stress  and  strain  in  the  test  article  will  become  nonlinear  when 
the  strain  corresponding  to  the  proportional  limit  is  exceeded.  In  the  working  models 
examined  herein  the  material  is  assumed  to  be  perfectly  elastic,  independently  of  the 
magnitude  of  strain. 

3.  Differences  in  constraint  conditions.  Rigid  end  fixtures  were  attached  to  the  free  end  of 
the  cylinder  and  the  nozzle,  as  described  in  Section  4.1.1.  Details  on  the  end  fixtures 
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are  not  given  in  the  ORNL  reports,  however  the  investigators  assumed  that  the  end 
fixtures  were  sufficiently  rigid  to  constrain  the  ends  of  the  cylinder  and  nozzle  so  as 
to  maintain  the  ends  as  plane  circles  [14] .  In  the  working  model  constructed  by  the 
ORNL  investigators  the  end  fixtures  were  represented  by  end  plates.  In  the  present 
investigation  the  fixture  attached  to  the  nozzle  was  represented  by  an  end  plate.  At  the 
free  end  of  the  cylinder  the  radial  and  tangential  displacement  components  were  set  to 
zero. 

4.  Differences  in  loading  conditions.  The  test  article  was  loaded  through  hydraulic  rams 
acting  on  the  end  fixtures.  The  accuracy  of  the  applied  load,  and  hence  the  accuracy 
of  the  stress  resultants,  was  determined  by  the  accuracy  of  the  load  cells.  The  transfer 
of  the  load  through  the  end  fixtures  was  through  mechanical  contact.  The  precise 
distribution  of  the  tractions  acting  on  the  ends  is  not  known.  In  the  working  models 
used  in  the  present  investigation  uniform  tractions  were  applied  on  the  free  end  of  the 
cylinder. 

5.  Symmetry.  The  working  model  was  assumed  to  be  perfectly  symmetric.  The  test  article 
was  not  perfectly  symmetric  and,  very  likely,  the  strain  gauges  were  not  installed  to  be 
perfectly  symmetric. 

It  is  seen  that  even  under  very  carefully  controlled  experimental  conditions  some  degree  of 
uncertainty  concerning  the  physical  system  is  present.  Some  of  these  uncertainties  can  be 
reduced,  other  uncertainties  either  cannot  be  reduced  or  may  not  be  feasible  to  reduce.  For 
example,  the  mean  value  of  the  elastic  constants  can  be  determined  by  simple  coupon  tests. 
The  dimensions  of  the  test  article  can  be  measured  with  high  accuracy.  On  the  other  hand,  it 
would  be  very  difficult  to  determine  the  distribution  of  the  tractions  or  constraint  conditions 
imposed  by  the  end  fixtures.  In  addition,  some  degree  of  uncertainty  is  associated  with  the 
instruments  employed  in  making  the  observations  and  the  effects  of  the  environment  on  the 
instruments. 

In  view  of  these  uncertainties  one  cannot  expect  very  close  correlation  between  computed 
and  experimental  data.  In  the  ORNL  experiments  the  largest  uncertainties  are  caused  by 
the  difficulties  associated  with  manufacturing  thin-walled  objects  to  tight  tolerances  and 
mathematical  representation  of  the  constraint  conditions. 

4.5  Numerical  results 

It  is  necessary  to  ascertain  that  the  errors  in  the  computed  data  are  well  below  the  threshold 
set  for  rejection  of  a  working  model.  In  this  investigation  the  following  steps  were  taken: 
(a)  The  relative  error  in  energy  norm  was  estimated.  This  provides  an  overall  view  of  the 
quality  of  the  approximation.  The  error  in  energy  norm  is  roughly  equivalent  to  the  RMS 
error  in  stresses  [7].  For  the  188-element  mesh  the  estimated  relative  error  in  energy  norm 
ranged  between  1.15  and  3.20  percent,  (b)  The  dependence  of  the  data  of  interest  on  the 
mesh  and  the  polynomial  degree  of  elements  was  examined.  It  was  found  that  the  data  of 
interest  are  substantially  independent  of  the  mesh  and  the  polynomial  degree  of  elements. 
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(c)  The  data  of  interest  were  examined  for  jump  discontinuities  at  inter-element  boundaries. 
Substantial  jump  discontinuities  in  locations  where  the  data  should  be  continuous  is  an 
indication  that  the  discretization  is  inadequate.  No  significant  discontinuities  were  found. 
These  are  necessary  conditions  that  an  approximate  solution  will  satisfy  when  the  errors  of 
approximation  are  not  large. 

Due  to  space  limitations,  only  some  of  these  procedures  are  illustrated  in  the  following. 
Details  of  the  finite  element  mesh  in  the  shell  intersection  region  are  shown  in  Fig.  13. 
The  size  of  the  shell  intersection  region  is  characterized  by  the  dimensions  ds  and  dn.  Unless 
otherwise  stated,  in  the  numerical  investigation  described  herein  these  dimensions  were  fixed: 
ds  =  0.25  in  (6.35  mm);  dn  =  0.20  in  (5.08  mm). 

The  layout  of  the  mesh  in  the  intersection  region  is  typical  of  meshes  used  in  ^extensions 
when  comer  singularities  are  present  [7].  Nodal  point  B  lies  on  the  line  of  intersection  between 
the  outer  surface  of  the  nozzle  and  the  outer  surface  of  the  shell  in  the  plane  of  symmetry. 
Nodal  point  A  lies  on  the  line  of  intersection  between  the  outer  surface  of  the  shell  and  the 
plane  of  symmetry. 

The  convergence  of  the  equivalent  stress  at  Point  A,  with  respect  to  the  number  of  degrees 
of  freedom  N,  is  shown  in  Fig.  14  on  a  semi-log  scale  for  the  fully  three-dimensional  model. 
It  is  seen  that  the  stress  is  substantially  independent  of  the  polynomial  degree  p  for  p  >  5. 
The  equivalent  stress  computed  from  the  finite  element  solution  cannot  be  close  to  its  exact 
value  if  this  criterion  is  not  satisfied. 


Figure  14:  Convergence  of  the  equivalent  stress  at  Node  A  shown  in  Fig.  13.  188-element  mesh. 

A  similar  plot  shown  in  Fig.  15  for  the  equivalent  stress  computed  for  point  B  clearly 
indicates  divergence.  This  is  caused  by  the  presence  of  an  edge  singularity.  The  equivalent 
stress  computed  for  the  exact  solution  of  the  elasticity  problem  is  infinity  in  this  point,  hence 
the  equivalent  stress  computed  from  the  finite  element  solution  cannot  converge  to  a  finite 
value.  Consequently,  extrapolation  of  strain  data  from  the  gauge  locations  to  the  intersection 
is  not  permissible. 
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Figure  15:  Convergence  of  the  equivalent  stress  at  Node  B  shown  in  Fig.  13.  188-element  mesh. 


Equivalent  stresses  in  strain  gauge  locations  nearest  to  the  intersection  region  in  the  plane 
of  symmetry  (toward  the  fixed-end  of  the  cylinder)  are  shown  in  Table  1.  These  data  were 
computed  from  the  solutions  obtained  for  the  fully  three-dimensional  model  with  p  ranging 
from  1  to  8.  The  number  of  degrees  of  freedom  (N)  are  shown  in  the  second  column.  Points 
C,  D,  E  and  F  are  shown  in  Fig.  13.  Each  stress  value  converges  strongly  to  a  limit  value, 
similar  to  the  convergence  shown  in  Fig.  14.  The  corresponding  experimental  data  and  the 
relative  differences  (DIFF),  using  the  computed  values  as  the  base,  are  shown  in  the  last  two 
rows.  It  is  seen  that  the  differences  are  very  substantial  in  three  of  the  four  points. 


Table  1:  EQUIVALENT  STRESSES  (PSI)  IN  THE  GAUGE  LOCATIONS  C,  D,  E,  F  IDENTIFIED  IN 
Fig.  13.  FULLY  THREE-DIMENSIONAL  MODEL.  THE  POINTS  ARE  LOCATED  ON  THE  FIXED-END 
SIDE  OF  THE  CYLINDER  (1  PSI=6.895  kPa). 


p 

N 

Pt.  C 

Pt.  D 

Pt.  E 

Pt.  F 

3 

6095 

11225 

11959 

13729 

16589 

4 

10828 

12461 

11680 

16729 

19687 

5 

17774, 

12455 

12158 

17644 

19441 

6 

27497 

12535 

12089 

17237 

19005 

7 

40561 

12542 

12021 

17108 

19104 

8 

57530 

12544 

11982 

17094 

18903 

EXPT 

15569 

14554 

16981 

13100 

DIFF.  (%) 

28.6 

24.1 

-0.3 

-28.9 

These  differences  are  attributed  primarily  to  variations  in  wall  thickness.  All  other  uncer¬ 
tainties  would  have  a  much  smaller  effect.  Inspection  of  the  test  article  revealed  variations  in 
wall  thickness  as  large  as  15  %  [14].  Variations  in  wall  thickness  affect  both  the  distribution 
and  magnitude  of  stresses.  The  magnitude  of  these  effects  depends  on  whether  the  solution 
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is  bending-  or  membrane-dominated.  Assuming  that  the  bending  moments  and  membrane 
forces  are  independent  of  the  wall  thickness,  in  a  bending  (resp.  membrane)  dominated  region 
10%  change  in  wall  thickness  causes  approximately  24%  (resp.  11%)  change  in  stress. 

In  experiments  designed  for  purposes  of  validation  it  is  necessary  to  eliminate  uncertain¬ 
ties  with  respect  to  the  object  of  the  experiment  as  much  as  possible  [30].  In  the  ORNL  test 
article  described  in  this  paper  the  dominant  source  of  uncertainty  is  the  unknown  variation 
in  wall  thickness.  Ideally,  the  external  and  internal  surfaces  would  be  carefully  measured 
and  the  actual  surfaces  would  be  used  in  the  working  model.  With  today’s  technology  it 
would  be  possible  to  produce  a  CAD  model  of  the  test  article  as  manufactured,  to  within 
0.001  in  (0.025  mm)  tolerance.  The  complexity  of  the  resulting  geometrical  description  of  the 
test  article  would  make  construction  of  the  working  model  significantly  more  complicated, 
however. 

In  the  immediate  vicinity  of  the  junction  there  are  large  bending  moments  that  decay 
with  respect  to  distance  from  the  junction.  The  computed  and  experimentally  obtained 
values  of  the  equivalent  stress  in  gauge  locations  along  the  line  of  intersection  of  the  plane  of 
symmetry  with  the  outer  surface  of  the  cylinder  on  the  fixed  end  side  are  tabulated  in  Table 
2.  It  is  seen  that  the  errors  are  larger  near  the  junction  than  away  from  it,  and  the  errors 
are  roughly  consistent  with  approximately  10%  variation  in  wall  thickness.  The  magnitude 
of  the  actual  variation  is  unknown. 

Table  2:  COMPUTED  (FEA)  AND  EXPERIMENTALLY  (EXP)  OBTAINED  VALUES  OF  THE  EQUIV¬ 
ALENT  STRESS  (PSI)  IN  GAUGE  LOCATIONS  ON  THE  INTERSECTION  OF  THE  PLANE  OF  SYM¬ 
METRY  WITH  THE  OUTSIDE  SURFACE  OF  THE  CYLINDER  AS  A  FUNCTION  OF  THE  DISTANCE 
(s)  FROM  THE  MID-SURFACE  OF  THE  NOZZLE  (INCHES). 


s 

FEA 

EXP 

err  (%) 

0.125 

12541 

15569 

24.14 

0.250 

10129 

12381 

22.23 

0.375 

8190 

8870 

8.31 

0.500 

6531 

7453 

14.13 

0.625 

5086 

5783 

13.70 

1.000 

2105 

2242 

6.49 

1.500 

1836 

1844 

0.42 

2.000 

2343 

2230 

-4.84 

3.000 

2431 

2292 

-5.71 

The  other  sources  of  error  are  inelastic  deformation  in  the  vicinity  of  the  junction,  vari¬ 
ations  in  material  properties,  errors  in  the  location  of  the  strain  rosettes,  errors  in  measure¬ 
ment,  errors  caused  by  idealization  of  the  actual  constraint  condition  and  errors  in  loading. 
None  of  these  errors  are  large  enough  to  explain  the  observed  differences.  For  example,  to 
test  the  effect  of  inelastic  deformation,  elastic  perfectly  plastic  material  behavior  was  as¬ 
sumed,  the  yield  stress  being  36.0  ksi  (248  MPa).  Using  the  deformation  theory  of  plasticity 
it  was  found  that  the  plastic  zone  was  confined  to  a  very  small  neighborhood  of  the  line  of 
intersection  of  the  outer  surfaces  and  hence  its  effect  on  the  stresses  in  the  gauge  locations 
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was  negligible.  For  a  discussion  on  the  other  effects  we  refer  to  [18]. 

A  comparison  between  data  computed  from  experimental  measurements  and  data  com¬ 
puted  from  the  fully  three-dimensional  working  model  along  the  inner  surface  of  the  cylinder, 
measured  from  the  mid-surface  of  the  nozzle,  is  shown  in  Fig.  16. 


Distance  from  the  intersection  [mm] 


Figure  16:  Equivalent  stress  on  the  inside  surface  of  the  cylinder  vs.  distance  from  the  mid-surface  of  the 
nozzle. 

Equivalent  stresses  computed  in  the  strain  gauge  locations  labelled  C,  D,  E,  F  in  Fig.  13 
(see  also  Table  1)  from  solutions  obtained  by  means  of  the  hierarchic  thin  solid  models  based 
on  the  spaces  S™9(Q, are  shown  in  Table  3.  The  estimated  relative  errors  in  energy  norm, 
denoted  by  (er)j g,  are  also  given  in  Table  3.  These  estimates  are  based  on  p-extension. 
Details  are  available  in  [7].  It  is  seen  that  the  stress  data  are  insensitive  to  the  choice  of  thin 
solid  model  characterized  by  q. 

Table  3:  EQUIVALENT  STRESSES  (PSI)  IN  THE  GAUGE  LOCATIONS  C,  D,  E,  F  IDENTIFIED  IN 
Fig.  13.  THIN  SOLID  MODELS.  THE  POINTS  ARE  LOCATED  ON  THE  FIXED-END  SIDE  OF  THE 
CYLINDER  (1  PSI=6.895  kPa). 


ds 

=  0.25  in  (6.35  mm),  dn 

=  0.2  in  (5.08  mm) 

p 

q 

(er)E 

N 

Pt.  C 

Pt.  D 

Pt.  E 

Pt.  F 

8 

l 

3.20 

41032 

12110 

11730 

16926 

18432 

8 

2 

1.15 

44766 

12570 

11978 

17094 

18904 

8 

3 

1.35 

47832 

12552 

11979 

17094 

18904 

8 

8 

1.35 

57530 

12544 

11982 

17094 

18903 

One  of  the  modeling  assumptions  is  the  selection  of  the  size  of  the  intersection  region 
indicated  in  Fig.  13  by  labels  ds  and  dn.  The  data  presented  so  far  was  computed  with 
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Table  4:  EQUIVALENT  STRESSES  (PSI)  IN  THE  GAUGE  LOCATIONS  C,  D,  E,  F  IDENTIFIED  IN 
Fig.  13.  THIN  SOLID  MODELS.  THE  POINTS  ARE  LOCATED  ON  THE  FIXED-END  SIDE  OF  THE 
CYLINDER  (1  PSI=6.895  kPa). 


ds 

=  0.13  in  (3.30  mm),  dn 

=  0.15  in  (3.81  mm) 

p 

q 

(er)E 

N 

Pt.  C 

Pt.  D 

Pt.  E 

Pt.  F 

8 

l 

2.44 

41032 

12111 

12470 

17076 

18169 

8 

2 

1.11 

44766 

12621 

12026 

17104 

19141 

8 

3 

1.32 

47832 

12564 

11941 

17104 

19143 

8 

8 

1.31 

57530 

12546 

11964 

17106 

19143 

ds  —  0.25  in  (6.35  mm)  dn  =  0.20  in  (5.08  mm).  Therefore  points  C,  D,  E,  F  were  located 
within  the  intersection  region.  Letting  ds  =  0.13  in  (3.30  mm)  dn  —  0.15  in  (3.81  mm)  the 
points  are  located  in  the  thin  solid  region.  Repeating  the  computations  the  results  shown  in 
Table  4  are  obtained.  On  comparing  Table  '3  with  Table  4  it  is  seen  that  the  differences  are 
negligibly  small.  This  indicates  that  the  working  models  are  insensitive  to  the  size  of  the 
intersection  region. 


4.6  Summary  and  conclusions 

The  intent  of  the  ORNL  investigation  was  to  determine  whether  the  finite  element  space 
characterized  by  the  mesh  shown  in  Fig.  11  and  the  HCT  plate  elements,  combined  with 
the  CLST  plane  stress  elements,  are  capable  of  approximating  strains  in  the  given  gauge 
locations.  We  understand  this  problem  to  consist  of  three  parts:  (a)  Whether  the  Novozhilov- 
Koiter  model  is  capable  of  representing  the  deformation  of  the  test  article  in  the  intersection 
region;  (b)  whether  the  three-dimensional  assembly  of  planar  plate  and  membrane  elements 
converge  to  the  exact  solution  of  the  Novozhilov-Koiter  model  as  the  size  of  the  elements  is 
reduced,  and  (c)  whether  the  mesh  shown  in  Fig.  11  is  suitable  for  controlling  the  errors  of 
discretization.  The  ORNL  investigation  did  not  address  these  questions  separately. 

In  the  present  investigation  the  need  to  address  these  questions  was  avoided  by  (a) 
demonstrating  that  hierarchic  thin  solid  formulations  provide  consistent  results  in  the  gauge 
locations  even  for  the  lowest  member  of  the  hierarchy,  characterized  by  q  =  1;  (b)  curved 
elements  mapped  by  smooth  functions  rather  than  planar  elements  were  used,  and  (c)  the 
data  of  interest  were  shown  to  be  substantially  independent  of  the  discretization.  Another 
important  difference  between  the  present  investigation  and  the  ORNL  investigation  is  that  in 
the  present  investigation  the  intersection  region  was  exempt  from  the  kinematic  assumptions 
of  the  thin  solid  formulation. 

The  ORNL  investigation  highlights  some  of  the  difficulties  and  limitations  of  experimen¬ 
tal  validation  of  working  models  for  shell  problems.  The  experimental  data  are  dominated 
by  uncertainties  caused  by  difficulties  associated  with  the  fabrication  of  thin- walled  objects 
to  exacting  tolerances.  Increasing  the  wall  thickness  would  reduce  errors  caused  by  manu¬ 
facturing  tolerances  but  then  the  thin  shell  model  may  not  be  applicable. 
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The  primary  objective  in  defining  a  working  model  is  to  account  for  all  physical  laws  and 
relationships  that  have  a  significant  influence  on  the  data  of  interest.  Therefore  the  choice 
of  a  working  model  depends  on  the  data  of  interest  and  has  to  be  validated  with  respect  to 
the  data  of  interest.  The  secondary  objective  is  to  identify  the  simplest  working  model  that 
will  satisfy  the  tolerances  set  for  the  data  of  interest. 

In  correlation  with  experimental  observations  a  working  model  is  tested  against  data 
that  are  either  observable  or  can  be  computed  from  observable  data.  In  the  model  problem 
discussed  herein  the  strain  components  in  surface  points  located  in  the  vicinity  of  the  shell 
intersection  were  measured  and  the  von  Mises  stresses  were  computed  in  the  gauge  locations. 

The  data  of  interest  are  generally  not  observable  and  cannot  be  computed  from  observable 
data.  For  example,  one  may  be  interested  in  the  maximum  value  of  the  integral  of  the  normal 
stress  over  some  small  area.  Therefore,  even  if  a  working  model  were  shown  to  be  successful  in 
predicting  certain  measured  data,  it  may  not  be  suitable  for  computing  other  data  of  interest. 
It  is  necessary  to  have  means  for  systematic  evaluation  of  the  effects  of  various  assumptions, 
incorporated  in  a  working  model,  on  the  data  of  interest.  Virtual  experimentation  based 
on  a  hierarchic  framework  of  models  and  hierarchic  discretizations  is  indispensable  in  the 
development  of  expert  knowledge. 

The  validity  of  a  working  model  cannot  be  established  by  experimental  correlation  in 
general.  The  purpose  of  validation  experiments  is  to  determine  whether  certain  necessary 
conditions  are  met  by  a  working  model.  Validation  experiments  cannot  establish  sufficient 
conditions  for  acceptance  of  a  working  model. 
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5  Generalized  plane  strain 


Dimensional  reduction  is  frequently  Used  in  engineering  practice.  The  goal  is  to  construct 
a  simplified  mathematical  model,  called  working  model,  so  that  the  exact  values  of  the 
data  of  interest  corresponding  to  the  simplified  model  are  substantially  the  same  as  those 
corresponding  to  the  fully  three-dimensional  model  it  is  expected  to  represent.  We  consider 
a  prismatic  body  of  length  £.  The  material  points  occupy  the  domain  D/,  defined  as  follows: 

Qe  =  {(x,  V,  z)  |  ( x ,  y )  €  u,  -i/2  <z<  i/2,  l  >  0}  (18) 

where  w  E  R2  is  a  bounded  domain  with  Lipschitz  boundary.  The  lateral  boundary  of  the 
body  is  denoted  by 


r*  =  {(x,  y,  z )  I  ( x ,  y)  e  du,  -i/2  <z<  i/2,  i  >  0}  (19) 

and  the  faces  are  denoted  by 

7±  =  {fo  lb  z)  |  (x,  y)eu,  z  =  ±i/2}.  (20) 

The  notation  is  shown  in  Fig.  17.  The  diameter  of  u  will  be  denoted  by  dw. 


Figure  17:  Notation. 

The  material  properties,  the  volume  forces  and  temperature  change  acting  on  €lt  and 
tractions  acting  on  T/  will  be  assumed  to  be  independent  of  z.  Therefore  the  x  y  plane  is 
a  plane  of  symmetry.  It  is  assumed  that  tractions  are  specified  on  the  entire  boundary  IV 
This  is  usually  called  the  first  fundamental  boundary  value  problem  of  elasticity. 

When  the  tractions  acting  on  7±  are  zero  and  i/du  «  1  then  such  problems  are  usually 
formulated  on  u  as  plane  stress  problems.  When  the  normal  displacements  and  shearing 
tractions  acting  on  7±  are  zero  then  such  problems  are  formulated  on  u  as  plane  strain 
problems,  independently  of  the  i/dw  ratio.  When  1  «  i/d^  and  the  tractions  acting  on  j± 
are  zero  then  the  three-dimensional  thermoelasticity  problem  on  is  called  the  generalized 
plane  strain  problem. 
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The  following  questions  were  addressed:  (a)  Under  what  conditions  is  the  plane  strain 
solution  asymptotically  correct;  (b)  what  corrections  to  the  plane  strain  solution,  that  can  be 
determined  by  solving  two-dimensional  problems,  are  necessary  to  obtain  an  asymptotically 
correct  solution  to  the  generalized  plane  strain  problem,  and  (c)  under  what  conditions  is 
the  plane  stress  solution  equivalent  to  the  generalized  plane  strain  solution. 

The  conditions  under  which  the  generalized  plane  strain  problem  is  asymptotically  (i.e., 
with  respect  to  £/du  — >  oo)  the  plane  strain  solution  were  stated  in  the  form  of  theorems.  It 
was  shown  also  that  in  special  cases  the  generalized  plane  strain  solution  is  the  plane  stress 
solution.  Details  and  examples  are  available  in  [31]. 
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